Fibroblasts as an in vitro model of circadian genetic and genomic studies

doi:10.21203/rs.3.rs-4272425/v1

Download PDF

Research Article

Fibroblasts as an in vitro model of circadian genetic and genomic studies

https://doi.org/10.21203/rs.3.rs-4272425/v1

This work is licensed under a CC BY 4.0 License

Journal Publication

published 01 Oct, 2022

Read the published version in European Neuropsychopharmacology →

Version 1

posted

You are reading this latest preprint version

Background Bipolar disorder (BD) is a heritable disorder characterized by shifts in mood that manifest in manic or depressive episodes. Clinical studies have identified abnormalities of the circadian system in BD patients as a hallmark of underlying pathophysiology. Fibroblasts are a well-established in vitro model for measuring circadian patterns. We set out to examine the underlying genetic architecture of circadian rhythm in fibroblasts, with the goal to assess its contribution to the polygenic nature of BD disease risk. We collected, from primary cell lines of 6 healthy individuals, temporal genomic features over a 48 hour period from transcriptomic data (RNA-seq) and open chromatin data (ATAC-seq).

Results The RNA-seq data showed that only a limited number of genes, primarily the known core clock genes such as ARNTL, CRY1, PER3, NR1D2 and TEF display circadian patterns of expression consistently across cell cultures. The ATAC-seq data identified that distinct transcription factor families, like those with the basic helix-loop-helix motif, were associated with regions that were increasing in accessibility over time. Whereas known glucocorticoid receptor target motifs were identified in those regions that were decreasing in accessibility.

Further evaluation of these regions using stratified linkage disequilibrium score regression analysis failed to identify a significant presence of them in the known genetic architecture of BD, and other psychiatric disorders or neurobehavioral traits in which the circadian rhythm is affected.

Conclusions In this study, we characterize the biological pathways that are activated in this in vitro circadian model, evaluating the relevance of these processes in the context of the genetic architecture of BD and other disorders, highlighting its limitations and future applications for circadian genomic studies.

Circadian rhythms

Transcriptomics

Epigenomics

Genomic Analysis

Neurobehavioral traits

It is estimated that the lifetime worldwide prevalence of bipolar disorder (BD) is 1% [1], with an estimated heritability of 60–85% [2]; [3]. Genome-wide association studies (GWAS) of BD are showing a highly polygenic genetic architecture of disease susceptibility with common genetic variants explaining 20% of the heritability [4]; [5]. BD is primarily characterized by shifts in mood, which result in manic or depressive episodes. Clinical studies have associated abnormalities of the circadian system in Bipolar disorder type 1 (BD1) patients as a hallmark component of its pathophysiology, with disturbed sleep quality being identified as an early symptom of manic episodes [6]. Furthermore, dysregulation of sleep and wake cycles during manic episodes include sleep abnormalities such as decrease in total sleep time, delta sleep, and REM latency [7]. These abnormalities have also extended to other circadian regulated systems such as cortisol levels. Both differences at morning levels of cortisol within BD subjects when compared to controls [8], as well as higher cortisol levels prior to a manic episode [9] have been reported. Despite these findings, the precise mechanisms of altered circadian rhythms in BD remain unclear.

The circadian rhythms synchronize physiological processes with the environment, creating and maintaining an internal 24 hour cycle. The main controller of the circadian cycle in mammals is the suprachiasmatic nucleus (SCN), a brain region located in the basal hypothalamus. It receives environmental cues, also called zeitgebers, such as light information from the retina which is relayed using synaptic and hormonal signaling [10] to the rest of the central nervous and peripheral systems. At the molecular level, the circadian machinery within every cell [11] consists of multiple transcriptional feedback loops, where core circadian genes BMAL1 and CLOCK induce the expression of their own repressors, PER1, PER2, PER3 and CRY1, CRY2. These genes modulate different layers of gene expression, from modifying the chromatin landscape to make certain regions of the genome more or less accessible [12], to post-transcriptional modifications altering the function of the associated proteins at specific times during the day [13]. Although disruptions in the circadian rhythms have been associated with neuropsychiatric traits, specifically in mood disorders [14], the direct interactions between them, as well as the contributions from genomic loci, are to be elucidated.

The localization of the SCN makes direct interaction and collection in humans impossible, with researchers instead using peripheral fibroblast cells to study the molecular and genetic components of this system [15]. These cells receive cortisol as a circadian signal from the SCN, through the hypothalamic-pituitary-adrenal axis (HPA). In order to study circadian rhythms using cell cultures, the cells need to be synchronized. One approach for this is treating the cells with dexamethasone, which elicits rhythm synchronization between the cells in a culture [15]. Dexamethasone binds to the glucocorticoid receptor, acting on the same pathways through which cortisol regulates circadian rhythms in vivo. This synchronization method has been used in conjunction with luciferase bioluminescence reporter assays to study the molecular dynamics of selected circadian genes in vitro [16]. Studies using these systems have been applied to both sleep disorders and BD. Although researchers were able to find differences in the period of expression of circadian genes in sleep disorders [17], similar studies using cells derived from BD1 patients were unable to detect significant [18] or replicable[19] differences.

Here we examine the genomic components of circadian related genetic regulation and general biology that this in vitro fibroblast model captures, and assess whether these features relate to the genetic architecture of BD susceptibility. For this, we collected longitudinal temporal sequencing data of both gene expression and accessible chromatin regions. The temporal gene expression was used to identify genes that display circadian oscillations and are under glucocorticoid control, as well as genes with distinct temporal patterns representing other biological pathways. The temporal accessible chromatin data was used to identify regions of the genome and associated transcription factor motifs that are implicated in the temporal regulation of gene expression. Finally, we examined whether the genomic regions showing temporal transcriptomic and epigenomic circadian profiles in primary cultures of fibroblasts were enriched in genetic association signals of BD or other related psychiatric and sleep-related phenotypes.

Temporal RNA-seq captures genes with distinct longitudinal expression patterns

Outside of the subset of genes that compose the core circadian transcriptional feedback loop, most rhythmic genes are tissue specific [20]. Within fibroblasts, we aimed to identify the overall longitudinal patterns of all the genes that are temporally regulated and classify them based on their temporal features. For this purpose, we collected RNA-seq data every 4 hours for a 48-hour period, from cell cultures of 6 human primary fibroblasts that were derived from a skin biopsy of subjects with no psychiatric disorders. To select these subjects, we confirmed that their cell lines displayed measurable circadian oscillations via a bioluminescence assay (Supplementary Fig. 2). After quality control, the temporal RNA-seq dataset consisted of n = 11,004 genes. We used a cubic spline regression model to identify genes that had a significant effect of time in their expression ([21]; [22]; [23]). This approach identified n = 2,767 (~ 25%) genes with significant evidence (False discovery rate (FDR) < 0.05) for temporal changes of gene expression levels. To cluster these genes according to distinct temporal patterns we applied the Weighted Gene Co-Expression Network Analysis (WGCNA)[24], which identifies genes with highly correlated expression levels. WGCNA produced 11 modules with eigengene values that captured the principal time patterns present in the expression of these genes (i.e: temporal modules; Supplementary Fig. 3). Gene ontology (GO) analysis of WGCNA modules with MetaScape [25] highlighted specific cell processes associated with distinct temporal patterns among 11 of these modules. Figure 1 depicts the eigengene values of 4 temporal modules with significant enrichment of GO terms (FDR adjusted by Benjamini-Hochberg method). Genes in the turquoise module, which show a linear decrease in expression over time, had GO terms for supramolecular fiber organization (p = 1e-15) and mRNA splicing via spliceosome (p = 2.5e-13). In compasion, genes in the blue module, which show a linear increase in expression, had a GO term for cellular response to hormone stimulus (p = 1.3e-9). The genes in the black module, which show an increase in expression that plateaus by the 16 hour time point (28 hours after dexamethasone treatment), were enriched for chromatin organization(p = 1e-12) and transcription elongation by RNA polymerase II (p = 7.9e-8) GO terms. The genes in the brown module, which show an expression pattern opposite of the black module, had a GO term for intracellular protein transport (p = 1e-18). Lastly, the purple module, which has genes with a peak expression at the 12 hour time point (24 hours after dexamethasone treatment), had GO terms for cell division (p = 1e-67) and mitotic cell cycle (p = 1e-60). The complete results of GO analysis for all the WGCNA modules are available in Supplementary Table 1.

WGCNA results did not yield a module of co-expressed genes with eigengene values representative of oscillating 24 hrs cycles resembling a circadian rhythm, nor were circadian related functional enrichment of GO terms found in any of the modules. Next we focused on closely inspecting known circadian genes for skin fibroblasts, identified by a previous in vivo array-based gene-expression circadian study on human skin cells [26]. Out of the 1,439 circadian genes reported in that study, we identified 267 genes in our dataset with significant changes in expression over time (Fig. 2A). Using the circadian detection tools JTK Cycle [27], LS [28], ARSER [29], Metacycle [30] and RAIN [31], we aimed to detect significant oscillations within these putative circadian genes in the complete temporal RNA-seq dataset. Among these methods, only JTK and ARSER identified significant periodic expression patterns (after Benjamini-Hochberg correction of 0.05) for the circadian gene NR1D2, and further only ARSER identified significant oscillations for 73 genes. However, the predicted period differed between the methods. For example, JTK predicted a period of 27.6 hours for NR1D2, while ARSER predicted 24.7 hours (Supplementary table 2 and supplementary files). Therefore instead of using these circadian detection tools, we applied smoothing-splines mixed effect models using the R package "sme" [32] to model the temporal features of these circadian genes (Fig. 2B and Supplementary Fig. 4). These models showed that for some of these circadian genes, such as CRY2 and NFIL3, the circadian expression pattern is only present in some of the cell cultures, whereas for genes such as NR1D2 and TEF, the circadian pattern is ubiquitous across cell cultures from different individuals. The fitted models for these circadian genes were then used in time warping analysis to group genes with known expression dynamics (Fig. 2C and Supplementary Fig. 4). From these expression patterns, we corroborate that NR1D2 expression follows its inhibition effect on ARNTL and CRY1 [33]. Similarly, PER3 expression follows its inhibition effect with ARNTL. Despite observing similar expression patterns in PER2 and PER3, these were not consistent across individuals (Supplementary Fig. 4). While the expected inhibition relationship between CRY1 and ARNTL was not present (Supplementary Fig. 5), this pattern of expression was also reported in the circadian dataset that was used as reference [26].

Temporal open chromatin levels measured by ATAC-seq highlights potential regulatory regions and transcription factor binding sites

To identify regions of the genome associated with the regulation and downstream effects of circadian genes, we collected ATAC-seq data following the same temporal design as with the RNA-seq dataset. Quality control metrics such as fraction of reads in peaks and transcription starting site for these samples is available in the supplementary material. After removing a cell line that did not pass quality control, we merged all overlapping regions of open chromatin, also known as peaks, across samples and time points as described previously [34], to define a common set of ATAC-seq signals (n = 126,057). We then used cubic spline regression models to identify peaks that have a significant change in accessibility over time. This approach yielded n = 7,568 (6%) time significant peaks, which were functionally annotated using ChipSeeker [35], a software that annotates peaks with the nearest gene and genomic regions (Supplementary Fig. 8B). Peaks with significant changes in accessibility over time showed a similar genomic distribution as the full dataset (Supplemetary Fig. 8A). Following the approach for the RNA-seq data, we applied WGCNA to cluster peaks with similar temporal patterns of accessibility changes (Supplementary Fig. 7).

WGCNA identified 4 different modules for the temporal patterns of chromatin accessibility, however the main pattern that characterizes these modules is an overall increase or decrease in accessibility. One module captured all the regions that were decreasing in accessibility (Fig. 3A), comprising 4,435 peaks. The other 3 modules showed regions increasing in accessibility. Individual motif enrichment analysis conducted with HOMER (Heinz et al. 2010, [36], showed similar enrichment across these modules, therefore we combined them into a single cluster of regions increasing in accessibility, in total 3,133 peaks. Regions that were decreasing in accessibility over time (Fig. 3A) had motif sequences for Fos (p = 1e-1047), Fra1 (p = 1e-1041), ATF3 (p = 1e-1034), BATF (p = 1e-1002), Fra2 (p = 1e-996), AP-1 (p = 1e-976), Jun-AP1 (p = 1e-681), Bach2 (p = 1e-330) and JunB (p = 1e-1001). Most of these transcription factors are part of the AP-1 transcription complex. Regions that were increasing in accessibility over time (Fig. 3B) had motif sequences for BHLHA15 (p = 1e-201), TCF4 (p = 1e-180), NeuroG2 (p = 1e-160), Twist2 (p = 1e-160), Pitx1 (p = 1e-186), Atoh1 (p = 1e-162), Tcf21(p = 1e-147), Olig2 (p = 1e-130), ZBTB18 (p = 1e-139) and NeuroD1 (p = 1e-123). These are dimerizing transcription factors that have the basic helix-loop-helix protein structural motif. In both types of regions HOMER identified the binding sequence of the glucocorticoid response element (GRE), although the rank for the GRE motif in regions that were decreasing in accessibility was higher. For the known circadian transcription factors, HOMER identified significant enrichment of the binding sequences for BMAL1 (p = 1e-29), NPAS2 (p = 1e-9), CLOCK (p = 1e-10), particularly within regions that had increasing accessibility over time. For the regions with decreasing accessibility over time, HOMER identified enrichment of NFIL3 (p = 1e-11).

Stratified Linkage Disequilibrium Score Regression (sLDSC) analysis

Functional annotation of the ATAC-seq dataset showed that approximately one third of the peak regions identified are located in distal intergenic regions, with unknown functions. Furthermore, it also showed that these regions displaying transient changes in chromatin state are located across the entire genome. To examine whether these open chromatin regions highlighted in our study are enriched for genetic susceptibility of BD and other neuropsychiatric traits, we used sLDSC (stratified linkage disequilibrium analysis)[37] to calculate the partitioned heritability of these features. For this approach we used published Psychiatric Genomics Consortium (PGC) summary statistics for BD [5], ADHD (Attention-Deficit/HyperactivityDisorder) [38], schizophrenia [39], PTSD (Post-traumatic stress disorder) [40], MDD (Major depression disorder) [41], insomnia [42], and the circadian trait of morningness [43]. We used the temporally significant ATAC-seq regions with 1 kilobases (kb) and 10 kb genomic windows in both downstream and upstream directions for each region. These ATAC-seq defined annotations were tested jointly with the baseline annotations included with sLDSC [37]. Figure 4 shows the enrichment for the traits tested from the ATAC-seq regions annotations as well as the baseline annotations (Full enrichment results are provided in the Supplementary Material). Among these, only the ATAC-seq regions that were decreasing in accessibility had a nominally significant (p value = 0.00463) less than expected presence for ADHD, and this effect was not present when the regions are extended by either 1kb or 10kb. In comparison, baseline annotations such as conserved regions in mammals showed a significant enrichment for all the traits (ADHD p value = 7.88e-11, schizophrenia p value = 1.67e-23, BD p value = 1.14e-8, MDD p value = 5.92e-22, insomnia p value = 2.65e-15, morningness p value = 3.05e-29), except PTSD (p value = 0.073). We did not identify significant enrichment of ATAC-seq regions in the other psychiatric and behavioral traits tested, indicating that these genomic regions with temporal trends in cromatin accessibility do not play a major role to their genetic architecture.

Cell cultures of peripheral tissues have been employed as models of in vitro circadian clock systems to study their molecular components [44] and the disorders in which they are disrupted [45]. For this particular model that uses primary cultures of fibroblasts derived from skin biopsies, we aimed to characterize the circadian features that are present at gene expression and chromatin accessibility levels, with the goal to identify the circadian genes that are engaged by this system as well as their associated regulatory genomic regions. From the longitudinal RNA-seq data we identified consistent circadian patterns of expression in a limited amount of genes such as ARNTL, CRY1, PER3, NR1D2 (Rev-erbBeta) and TEF, but observed noticeable differences in the expression patterns between cell cultures. When compared to a recent in vivo circadian human skin dataset [26], we identified 267 out of the 1,439 circadian genes previously identified in this tissue to have a significant effect of time in their expression. This limited overlap could indicate that this in vitro model for studying circadian rhythms is constrained to the circadian genes that are directly activated by a glucocorticoid-like stimulus. Glucocorticoid response elements have been identified for circadian genes such as PER1, PER2, PER3, CRY1, CRY2, Rev-erbAlpha (NR1D1), Rev-erbBeta (NR1D2), DBP, NPAS2 and BMAL1 [46]. Consistent with these results we identified circadian rhythmicity in most of those genes that had previously been identified to have a glucocorticoid response element (GRE), with the exception of NR1D1. While we do identify expression levels from NR1D1, the lack of a significant circadian oscillation in comparison to the strong results from NR1D2 could be consistent with their known redundant function for circadian rhythms [47].

We found robust circadian expression patterns for NR1D2. However, differences in the lenght of the predicted periods across tools indicates that estimating period duration from longitudinal RNA-seq data is not a straightforward problem. This could be due to the the small number of subjects used, leading to insufficient power. Furthermore, while ARSER identified 73 genes with significant periodic expression patterns, this software is known to have a high false positive rate in high resolution data [30]. Interestingly, while the expression pattern of CRY1 follows the expected inhibition by NR1D2 [48], it does not reflect the expected inhibitory action on ARNTL, nor the similar phase pattern with its heterodimer partner PER3. For the relationship with PER3, CRY1 has been previously reported to have a known phase delay with the PER genes (PER1,2,3)[49], which could be attributed to the multiple binding sites that CRY1 has for different circadian modulators, resulting in stimulus and tissue specific temporal dynamics. Based on the time patterns, changes in the expression of CRY1 appear to precede the expression of ARNTL by 4 to 8 hours. Similar expression patterns between these genes were reported in a previous in vivo study of human skin cells [26], indicating that this model was able to replicate some of the circadian temporal dynamics seen within tissue.

ATAC-seq data can provide an untargeted yet comprehensive view of chromatin accessibility changes over time. Within this fibroblast in vitro model, we mainly identified genomic regions with linear increases and decreases of chromatin accessibility. Although we did not found circadian temporal patterns in chromatin accessibility, previous studies conducted in vivo have reported such patterns [50]. Specifically, proteins such as CLOCK and BMAL1 have been found to associate and interact with chromatin remodeling and chromatin modifying enzymes [51], as well as act as pioneer factors by directly modifying chromatin accessibility [12]. The absence of anticipated oscillations in chromatin within our dataset, as opposed to our observations in genes, could be due to multiple reasons. One possibility is that the mechanisms governing oscillatory chromatin changes may be exclusive to in vivo conditions. Under a physiological setting, cells within a tissue are exposed to multiple stimuli that act as Zeitgebers, such as sunlight (exposure), metabolic signals, temperature, and hormones like cortisol [52]. The exposure to these signals is under a rhythmic control, with levels cycling throughout the day [53]. By using a single exposure to dexamethasone in this model, we are missing the cyclic aspect of cortisol response present in vivo, as well as other effects that could be due to the coupling of Zeitgebers. However, our ATAC-seq dataset does replicate previous findings on the broader role of glucocorticoids in the chromatin landscape. Within regions with decreasing accessibility post-dexamethasone, the motif enrichment analyses identified the motifs for the GRE as well as for members of the AP-1 transcriptional complex. These regions may have initially opened due to the dexamethasone treatment (for synchronization of the cells), but are closing without further continued exposure of dexamethasone. Regions with increasing accessibility may have initially closed due to the dexamethasone treatment, and this is consistent with motif profiles that were unrelated to direct glucocorticoid receptor (GR) targets [54]. The dataset, however, lacks a Chip-seq analysis for GR occupancy, and we are limited to confirm if the identified regions are indeed due to GR activity. The strong glucocorticoid effects observed in our data underscore the need for further exploration of circadian influences on chromatin regulation in fibroblast cell culture models. Other methods for synchronizing these cells and studying the circadian rhythms, such as a switch to serum free media [15], Forskolin treatment [55], or temperature cycles [56], could result in different types of chromatin regulation and gene expression dynamics. This study raises questions about the context-dependent nature of chromatin remodeling events and emphasizes the need to evaluate different synchronization methods to ascertain their implications for circadian rhythms.

The ATAC-seq data showed genome wide transient changes in chromatin conformation, with most of these changes occurring within regions of unknown functions. To evaluate the relevance of these genomic regions for BD and other psychiatric traits, we used partitioned sLDSC regression. This tool identifies genetic susceptibility enrichment for a particular trait across the whole genome and within specific genomic annotations. The partitioned sLDSC analysis mainly showed a significant deflation for the chromatin regions that were decreasing in accessibility over time with ADHD. Although not significant, it mirrored the enrichment for the regions that were increasing in accessibility. When expanding the genomic regions by either 1kb or 10kb both the magnitude and the significance of the enrichment are lost, indicating that this effect could be highly localized for these regions. For the other traits that we examined, the enrichment from the ATAC-seq regions were also attenuated when the genomic region was extended. These results show that these regions with linear changes in chromatin accessibility identified here may not play a relevant role for these neuropsychiatric traits. However the attenuation observed when expanding the genomic window suggests that any relevant signal may be specific to those genomic positions.

The lack of an overlap between the temporal regulatory regions identified in this study and the known genetic architecture of BD could indicate four different interpretations. First, although peripheral tissues capture the genome of an individual, they don’t recapitulate brain molecular physiology, the main tissue implicated in the pathophysiology of BD. Second, it could be that the disruptions in the circadian rhythm are not under strong genetic control and are actually influenced by other downstream processes, such as post-translational modifications and differences at the protein level. Third, the specific circadian pathways that are engaged in this in vitro model by dexamethasone (i.e: glucocorticoids) are not part of the genetic architecture of BD. This however does not exclude other circadian pathways that could be engaged by a different synchronizing stimulus, such as serum [57], [58]. Furthermore, there are other stimuli that also act in different ways with the circadian system, such as temperature. Whereas dexamethasone acts through binding of the glucocorticoid receptor, temperature affects heat-shock proteins [56]. Lastly, the disrupted circadian phenotype is an episodic state in BD patients, not a constant trait. This could indicate that rather than the regular circadian system being affected by BD, it is the ability to deal with circadian stressors and disruptors that is implicated in BD disease susceptibility.

The circadian analysis of gene expression and chromatin accessibility data faced limitations that could be attributed to variability among cell lines from different subjects. Human skin cell studies, both in vivo [26] and in vitro [59], have demonstrated that genetic differences contribute to variations in circadian gene expression's phase and amplitude. The variability observed in this dataset was therefore not unexpected, but remains a factor to be considered for this kind of studies. Another potential source of variability was the data collection scheme, involving 13 separate cell cultures for each individual cell line. Distinctions in cell cycle state and growth rates among these cultures might have influenced the data. Previous research has shown that cell cycle and circadian rhythms are coupled processes([60], [61]), and that these rhythms can be impacted by cell density [62]. Our approach, utilizing a 5% FBS culture that minimizes cell growth, aimed to control for both of these factors. Our lab's prior work also confirmed that the cell density used for our study allows for the production of rhythmic circadian cycles in these cells(Supplementary Fig. 2). Although various factors known to influence circadian rhythms could have contributed to the variability in this dataset, certain circadian genes appeared resilient, consistently producing rhythmic cycles across cultures and individual cell lines. This could highlight specific genes' resilience to various sources of variation in this kind of studies.

Cell lines and Culture

Fibroblasts were isolated by taking skin biopsies from the nether region from subjects without known psychiatric disorders. Fibroblast cultures were established following standard procedures [63] and stored as frozen aliquots in liquid nitrogen. 6 fibroblast cell lines matched for sex, age and passage number were thawed out and grown to confluence in T75 culture flasks in standard culture media (DMEM containing 10% fetal bovine serum (FBS) and 1x Penicillin-Streptomycin).

Upon reaching confluence, 5x10^4 cells were plated per line into 13 different 6 well plates (1 well per line per plate). All 6 lines were collected in the same experiment for the RNA-seq experiment. Due to the labor-intensive nature of the ATAC protocol and the need to process cells fresh, the 6 lines were split into 2 batches, so 3 lines per batch were processed.

Assessment of Circadian Expression in vitro

In order to collect RNA or cells every 4 hours for 48 hours, cells were split into two batches, which were reset 12 hours apart (see supplementary Fig. 9). Cells were reset 12 hours before the first collection to exclude the acute effects of dexamethasone and variation in synchronization conditions [59]. 5 days after being plated the cells from batch one were synchronized by treatment with 100 nM Dexamethasone for 30 min. Cells were then washed with PBS and switched to collection media (DMEM containing 5% FBS and 1x Penicillin-Streptomycin). Lower concentration of FBS was used in this media to stop the cells from growing during the experiment, in order to keep all time points at approximately the same culture density. 12 hours later cells from batch 2 were synchronized and switched to collection media and the RNA/cell collection was started (from batch one).

RNA and Cell collection

For the collection of RNA, cells were lysed using 350uL RLT lysis buffer from the Qiagen RNeasy mini kit. Lysed cells were then scraped off the plate, transferred to a Qiaschredder (Qiagen 79656) and centrifuged for 2 min at max speed to further homogenize. Cell lysates were kept in -80 until extraction.

For the collection of cells for the ATAC protocol, cells were dissociated using 500uL of prewarmed TrypLE (ThermoFisher 12604013) and left for 5 min at 37℃. TrypLE was inactivated using 500uL of DMEM. Cells were then counted using the Logos Biosystems LUNA-FL automated cell counter, and 50x10^4 cells were used as input for tagmentation. Tagmented DNA for library preparation was collected following the previously described protocol [64].

RNA extraction

RNA from cell lysates was extracted using the Qiagen RNeasy mini kit (Qiagen 74106). Cell lysates were extracted in a randomized order to prevent batch effects in downstream analysis. In order to collect total RNA including small RNAs, the standard extraction protocol (Purification of Total RNA from Animal Cells using Spin Technology) was adjusted by making the following changes: (i) adding 1.5 volumes of 100% ethanol, instead of 70%, after the lysis step (step 4 in handbook protocol) and (ii) adding 700 mL of buffer RWT (Qiagen 1067933) instead of the provided RW1 (step 6 in handbook protocol).

RNA and ATAC sequencing

For the RNA sequencing, library preps were made using the Lexogen QuantSeq 3’ mRNA-Seq Library Prep Kit and sequenced with 65-base single end reads, and sequenced at a targeted depth of 3.8M reads per sample, which is well above the recommended minimum 1M reads per sample read depth for these types of libraries. ATAC seq libraries were generated following the previously described protocol [64] and sequenced with 75-base double end reads, and sequenced at a targeted depth of 61M reads per sample. Library preparation and sequencing was performed at the UCLA Neuroscience Genomics Core (https://www.semel.ucla.edu/ungc). All samples were sequenced on a Illumina HiSeq 4000 sequencer.

RNA-seq data processing and analysis

Fastqc [65] software was used to assess the quality of the read files. Low quality reads were trimmed using TrimGalore and Cutadapt.

Alignment of reads was performed with the STAR[66] software and to human gene ensembl version GrCh38. STAR was indexed to the genome using the –runMode genomeGenerate function. For aligning, STAR was run with the parameters –outFilterType BySJout –outFilterMultimapNmax 20 –alignSJoverhangMin 8 –alignSJDBoverhangMin 1 –outFilterMismatchNmax 999 --outFilterMismatchNoverLmax 0.1 --alignIntronMin 20 --alignIntronMax 1000000 --alignMatesGapMax 1000000.

Samtools was used to index the aligned files from STAR.

Read counts were associated with genes using featureCounts software with the NCBI GRCh38 gene annotation file.

Analysis of the RNA-seq data used the R packages limma, Glimma and edgeR, as previously described (Law et al., 2018). Genes with low read counts were removed and reads were normalized by CPM. GeneIDs were converted to Gene Symbols using the package Homo.sapiens.

WGCNA [24] software was used to classify genes with similar temporal expression patterns. WGCNA was run using a power value of 12 obtained from diagnostic plots and with the "signed" argument. MetaScape was used for Gene Ontology analysis of the resulting gene sets from WGCNA.

Following the method described in [67], MetaCycle [30] was used to run circadian detection tools such as ARSER, JTK [27], LS and metacycle. In order to integrate results from the different individuals, the function meta3d was used from the Metacycle R package. RAIN [31] was run separately.

According to the EdgeR user guide, cubic splines were generated using the splines package in R, with the ns function and 5 degrees of freedom. Resulting p-values were corrected using a false discovery rate of 0.05. Significant genes were then compared to a previously published dataset of circadian human skin gene expression, resulting in 267 genes that were classified according to their time series using the “dtwclust” R package. The resulting clusters are available in the supplementary material. This analysis was also performed with only the known circadian clock genes that had consistent expression patterns across the cell-lines.

ATAC-seq data processing and analysis

Sequenced open chromatin data from the ATAC-seq assay followed the standard ENCODE Pipeline for the identification of open chromatin regions (OCRs) of the genome. The steps included using fastqc to evaluate the quality of the sequenced library. Followed by trimming of low quality reads with Trimgalore and Cutadapt. Alignment of the raw reads data to human gene Ensembl version GRCh38 was performed using bowtie2 [68] with a 2kb insert size and allowing up to 4 alignments. Reads within black-listed regions alongside PCR duplicates were removed with samtools. MACS2 [69] software was used to identify OCRs with parameters -g hs -q 0.01 –nomodel –shift − 100 –extsize 200 –keep-dup all -B. PCR. Quality control metrics for the ATAC-seq dataset such as peak counts, PCR bottlenecking coefficients, fraction of reads in peaks and enrichment of transcription starting site are provided in the supplementary material.

To compare the ATAC-seq signal across timepoints and subjects, we created a consensus bed file using the bedtools[70] function merge function, combining all the overlapping peak regions across timepoints and subjects into a single file. The Featurecounts software was then used to assign read counts to those regions. Read counts were normalized by RPM.

WGCNA [24] software was used to classify peaks with similar temporal accessibility patterns. WGCNA was run using a power value of 12 obtained from diagnostic plots and with the "signed" argument, identifying 4 modules after merging.

HOMER (Heinz et al. 2010) findMotifsGenome.pllp program was used to identify enriched transcription factor motifs first individually in the peaks that belonged to the largest WGCNA modules, as well as in the resulting set of grouping all the modules that displayed a similar increasing or decreasing pattern of accessibility.

Stratified linkage disequilibrium score regression analysis

Following the procedure described in [71], we applied an extension to sLDSR, a statistical method that partitions SNP-based heritability(h2) from GWAS summary statistics [37]. We ran sLDSR (ldsc.py –h2), using an ancestry-match 1000 Genomes Project phase 3 release reference panel, for each annotation of interest while accounting for the full baseline model, as recommended by the developers ([37]; [72]), and an extra annotation of all the ATAC-seq detected in our in vitro model (n = 3126 for peaks that were decreasing in accessibility, n = 4415 for peaks increasing in accessibility), as well as extension of these regions by 1kb and 10kb genomic windows in both directions.

With the knowledge of the specific features that this in vitro model is able to capture of the circadian system, we advise care when interpreting the results of such experiments, as they may be heavily influenced by genetics, cell culture factors, and, crucially, the circadian cycle synchronization method. This can, inadvertently, lead to a narrowing in the scope of studies of circadian rhythms in the context of neuropsychiatric traits. While the biology that this model captures after circadian synchronization induced by dexamethasone treatment does not seem to be directly involved in the known genetic architecture of BD, this model can still be applied to scientific questions that cannot be explored directly in human subjects. For instance, this model could be employed to characterize the specific biological pathways that are engaged during circadian distress. Notably, the dysregulated circadian phenotype in BD patients is characterized by episodic events rather than a static state, emphasizing the dynamic nature of the subjects. Responses to circadian distress could be directly compared between fibroblast cell lines derived from BD patients and healthy subjects. Additionally, the accessibility that this in vitro model provides could be used to study the effect of lithium, the most commonly used prescribed drug treatment for BD, during such circadian distress.

Ethics approval and consent to participate

Written informed consent was obtained from participants and approved by institutional review boards at the University of California Los Angeles

Consent for publication

Not applicable

Availability of data and materials

The datasets supporting the conclusions of this article are available in the https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE263711 repository and the https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE263713 repository.

Competing interests

The authors declare that they have no competing interests

Funding

This study was supported by funding from the National Institutes of Health (NIH), research grants R01 MH090553, R01 MH115676, the NARSAD Distinguished Investigator Grant, (to Roel Ophoff) and NS048004 T32 Training grant in Neurobehavioral genetics.

Author’s contributions

MF carried out the statistical analysis, interpretation and discussion of the results, and drafted the manuscript. MeB contributed to the conception, design and execution of the experimental protocol, grew cell cultures and collected the genomic data, and was involved in the revision of the manuscript. TB and JH helped with designing the statistical analysis, and substantively reviewed the manuscript. MaB and RK substantively reviewed the manuscript. RO contributed to the conception and design of experimental protocol, discussion of the results. All authors have reviewed and approved the manuscript.

Acknowledgments

We would like to acknowledge the advice and expertise of Dr. Christopher Colwell for this project, and to Dr. Annet van Bergen for their work in collecting the tissue samples used in this study. We thank the donors for their willingness to provide a skin biopsy for the generation of the primary cell cultures.

Moreira ALR, Van Meter A, Genzlinger J, Youngstrom EA. Review and Meta-Analysis of Epidemiologic Studies of Adult Bipolar Disorder. J Clin Psychiatry. 2017;78:e1259–69.
Song J, Bergen SE, Kuja-Halkola R, Larsson H, Landén M, Lichtenstein P. Bipolar disorder and its relation to major psychiatric disorders: a family-based study in the Swedish population. Bipolar Disord. 2015;17:184–93.
Bienvenu OJ, Davydow DS, Kendler KS. Psychiatric diseases versus behavioral disorders and degree of genetic influence. Psychol Med. 2011;41:33–40.
Stahl EA, Breen G, Forstner AJ, et al. Genome-wide association study identifies 30 loci associated with bipolar disorder. Nat Genet. 2019;51:793–803.
Mullins N, Forstner AJ, O’Connell KS, et al. Genome-wide association study of more than 40,000 bipolar disorder cases provides new insights into the underlying biology. Nat Genet. 2021;53:817–29.
Leibenluft E, Albert PS, Rosenthal NE, Wehr TA. Relationship between sleep and mood in patients with rapid-cycling bipolar disorder. Psychiatry Res. 1996;63:161–8.
Levenson J, Frank E. Sleep and circadian rhythm abnormalities in the pathophysiology of bipolar disorder. Curr Top Behav Neurosci. 2011;5:247–62.
Girshkin L, Matheson SL, Shepherd AM, Green MJ. Morning cortisol levels in schizophrenia and bipolar disorder: a meta-analysis. Psychoneuroendocrinology. 2014;49:187–206.
van den Berg MT, Wester VL, Vreeker A, Koenders MA, Boks MP, van Rossum EFC, Spijker AT. Higher cortisol levels may proceed a manic episode and are related to disease severity in patients with bipolar disorder. Psychoneuroendocrinology. 2020;119:104658.
Le Minh N, Damiola F, Tronche F, Schütz G, Schibler U. Glucocorticoid hormones inhibit food-induced phase-shifting of peripheral circadian oscillators. EMBO J. 2001;20:7128–36.
Schibler U, Sassone-Corsi P. A web of circadian pacemakers. Cell. 2002;111:919–22.
Menet JS, Pescatore S, Rosbash M. CLOCK:BMAL1 is a pioneer-like transcription factor. Genes Dev. 2014;28:8–13.
Robles MS, Humphrey SJ, Mann M. Phosphorylation Is a Central Mechanism for Circadian Control of Metabolism and Physiology. Cell Metab. 2017;25:118–27.
Walker WH 2nd, Walton JC, DeVries AC, Nelson RJ. Circadian rhythm disruption and mental health. Transl Psychiatry. 2020;10:28.
Yamazaki S, Takahashi JS. Real-time luminescence reporting of circadian gene expression in mammals. Methods Enzymol. 2005;393:288–301.
Nakahata Y, Akashi M, Trcka D, Yasuda A, Takumi T. The in vitro real-time oscillation monitoring system identifies potential entrainment factors for circadian clocks. BMC Mol Biol. 2006;7:5.
Hida A, Ohsawa Y, Kitamura S, et al. Evaluation of circadian phenotypes utilizing fibroblasts from patients with circadian rhythm sleep disorders. Transl Psychiatry. 2017;7:e1106.
Yang S, Van Dongen HPA, Wang K, Berrettini W, Bućan M. Assessment of circadian function in fibroblasts of patients with bipolar disorder. Mol Psychiatry. 2009;14:143–55.
McCarthy MJ, Wei H, Marnoy Z, Darvish RM, McPhie DL, Cohen BM, Welsh DK. Genetic and clinical factors predict lithium’s effects on PER2 gene expression rhythms in cells from bipolar disorder patients. Transl Psychiatry. 2013;3:e318.
Zhang R, Lahens NF, Ballance HI, Hughes ME, Hogenesch JB. A circadian gene expression atlas in mammals: implications for biology and medicine. Proc Natl Acad Sci U S A. 2014;111:16219–24.
Wang Y, Ke C, Brown MB. Shape-invariant modeling of circadian rhythms with random effects and smoothing spline ANOVA decompositions. Biometrics. 2003;59:804–12.
Qin L, Guo W. Functional mixed-effects model for periodic data. Biostatistics. 2006;7:225–34.
Madden JM, Li X, Kearney PM, Tilling K, Fitzgerald AP. Exploring diurnal variation using piecewise linear splines: an example using blood pressure. Emerg Themes Epidemiol. 2017;14:1.
Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9:559.
Zhou Y, Zhou B, Pache L, Chang M, Khodabakhshi AH, Tanaseichuk O, Benner C, Chanda SK. Metascape provides a biologist-oriented resource for the analysis of systems-level datasets. Nat Commun. 2019;10:1523.
Del Olmo M, Spörl F, Korge S, et al. Inter-layer and inter-subject variability of diurnal gene expression in human skin. NAR Genom Bioinform. 2022;4:lqac097.
Hughes ME, Hogenesch JB, Kornacker K. JTK_CYCLE: an efficient nonparametric algorithm for detecting rhythmic components in genome-scale data sets. J Biol Rhythms. 2010;25:372–80.
Glynn EF, Chen J, Mushegian AR. Detecting periodic patterns in unevenly spaced gene expression time series using Lomb-Scargle periodograms. Bioinformatics. 2006;22:310–6.
Yang R, Su Z. Analyzing circadian expression data by harmonic regression based on autoregressive spectral estimation. Bioinformatics. 2010;26:i168–74.
Wu G, Anafi RC, Hughes ME, Kornacker K, Hogenesch JB. MetaCycle: an integrated R package to evaluate periodicity in large scale data. Bioinformatics. 2016;32:3351–3.
Thaben PF, Westermark PO. Detecting rhythms in time series with RAIN. J Biol Rhythms. 2014;29:391–400.
Berk M. (2018) Smoothing-splines Mixed-effects Models in R using the sme Package: a Tutorial. https://rdrr.io/cran/sme/f/inst/doc/Tutorial.pdf. Accessed 10 Jun 2022.
Rijo-Ferreira F, Takahashi JS. Genomics of circadian rhythms in health and disease. Genome Med. 2019;11:82.
Keele GR, Quach BC, Israel JW, et al. Integrative QTL analysis of gene expression and chromatin accessibility identifies multi-tissue patterns of genetic regulation. PLoS Genet. 2020;16:e1008537.
Yu G, Wang L-G, He Q-Y. ChIPseeker: an R/Bioconductor package for ChIP peak annotation, comparison and visualization. Bioinformatics. 2015;31:2382–3.
Yan F, Powell DR, Curtis DJ, Wong NC. From reads to insight: a hitchhiker’s guide to ATAC-seq data analysis. Genome Biol. 2020;21:22.
Finucane HK, Bulik-Sullivan B, Gusev A, et al. Partitioning heritability by functional annotation using genome-wide association summary statistics. Nat Genet. 2015;47:1228–35.
Demontis D, Walters RK, Martin J, et al. Discovery of the first genome-wide significant risk loci for attention deficit/hyperactivity disorder. Nat Genet. 2019;51:63–75.
Trubetskoy V, Pardiñas AF, Qi T, et al. Mapping genomic loci implicates genes and synaptic biology in schizophrenia. Nature. 2022;604:502–8.
Nievergelt CM, Maihofer AX, Klengel T, et al. International meta-analysis of PTSD genome-wide association studies identifies sex- and ancestry-specific genetic risk loci. Nat Commun. 2019;10:4558.
Howard DM, Adams MJ, Clarke T-K, et al. Genome-wide meta-analysis of depression identifies 102 independent variants and highlights the importance of the prefrontal brain regions. Nat Neurosci. 2019;22:343–52.
Watanabe K, Jansen PR, Savage JE, et al. Genome-wide meta-analysis of insomnia prioritizes genes associated with metabolic and psychiatric pathways. Nat Genet. 2022;54:1125–32.
Jones SE, Tyrrell J, Wood AR, et al. Genome-Wide Association Analyses in 128,266 Individuals Identifies New Morningness and Sleep Duration Loci. PLoS Genet. 2016;12:e1006125.
Balsalobre A, Marcacci L, Schibler U. Multiple signaling pathways elicit circadian gene expression in cultured Rat-1 fibroblasts. Curr Biol. 2000;10:1291–4.
Kripke DF, Nievergelt CM, Joo E, Shekhtman T, Kelsoe JR. Circadian polymorphisms associated with affective disorders. J Circadian Rhythms. 2009;7:2.
So AY-L, Bernal TU, Pillsbury ML, Yamamoto KR, Feldman BJ. Glucocorticoid regulation of the circadian clock modulates glucose homeostasis. Proc Natl Acad Sci U S A. 2009;106:17582–7.
Liu AC, Tran HG, Zhang EE, Priest AA, Welsh DK, Kay SA. Redundant function of REV-ERBalpha and beta and non-essential role for Bmal1 cycling in transcriptional regulation of intracellular circadian rhythms. PLoS Genet. 2008;4:e1000023.
Chiou Y-Y, Yang Y, Rashid N, Ye R, Selby CP, Sancar A. Mammalian Period represses and de-represses transcription by displacing CLOCK-BMAL1 from promoters in a Cryptochrome-dependent manner. Proc Natl Acad Sci U S A. 2016;113:E6072–9.
Fustin JM, O’Neill JS, Hastings MH, Hazlerigg DG, Dardente H. Cry1 circadian phase in vitro: wrapped up with an E-box. J Biol Rhythms. 2009;24:16–24.
Koike N, Yoo S-H, Huang H-C, Kumar V, Lee C, Kim T-K, Takahashi JS. Transcriptional architecture and chromatin landscape of the core circadian clock in mammals. Science. 2012;338:349–54.
Zhu Q, Belden WJ. Molecular Regulation of Circadian Chromatin. J Mol Biol. 2020;432:3466–82.
Roenneberg T, Merrow M. The Circadian Clock and Human Health. Curr Biol. 2016;26:R432–43.
Chauhan S, Norbury R, Faßbender KC, Ettinger U, Kumari V. Beyond sleep: A multidimensional model of chronotype. Neurosci Biobehav Rev. 2023;148:105114.
Bothe M, Buschow R, Meijsing SH. Glucocorticoid signaling induces transcriptional memory and universally reversible chromatin changes. Life Sci Alliance. 2021. https://doi.org/10.26508/lsa.202101080.
Yagita K, Okamura H. Forskolin induces circadian gene expression of rPer1, rPer2 and dbp in mammalian rat-1 fibroblasts. FEBS Lett. 2000;465:79–82.
Saini C, Morf J, Stratmann M, Gos P, Schibler U. Simulated body temperature rhythms reveal the phase-shifting behavior and plasticity of mammalian circadian oscillators. Genes Dev. 2012;26:567–80.
Balsalobre A, Damiola F, Schibler U. A serum shock induces circadian gene expression in mammalian tissue culture cells. Cell. 1998;93:929–37.
Iyer VR, Eisen MB, Ross DT, et al. The Transcriptional Program in the Response of Human Fibroblasts to Serum. Science. 1999;283:83–7.
Brown SA, Fleury-Olela F, Nagoshi E, Hauser C, Juge C, Meier CA, Chicheportiche R, Dayer J-M, Albrecht U, Schibler U. The period length of fibroblast circadian gene expression varies widely among human individuals. PLoS Biol. 2005;3:e338.
Nagoshi E, Saini C, Bauer C, Laroche T, Naef F, Schibler U. Circadian gene expression in individual fibroblasts: cell-autonomous and self-sustained oscillators pass time to daughter cells. Cell. 2004;119:693–705.
Farshadi E, van der Horst GTJ, Chaves I. Molecular Links between the Circadian Clock and the Cell Cycle. J Mol Biol. 2020;432:3515–24.
Noguchi T, Wang LL, Welsh DK. Fibroblast PER2 circadian rhythmicity depends on cell density. J Biol Rhythms. 2013;28:183–92.
Villegas J, McPhaul M. (2005) Establishment and culture of human skin fibroblasts. Curr Protoc Mol Biol Chap. 28:Unit 28.3.
Buenrostro JD, Wu B, Chang HY, Greenleaf WJ. (2015) ATAC-seq: A Method for Assaying Chromatin Accessibility Genome-Wide. Curr Protoc Mol Biol 109:21.29.1–21.29.9.
Andrews S. (2010) FastQC A Quality Control tool for High Throughput Sequence Data. http://www.bioinformatics.babraham.ac.uk/projects/fastqc/. Accessed 1 Apr 2024.
Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, Batut P, Chaisson M, Gingeras TR. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29:15–21.
Mei W, Jiang Z, Chen Y, Chen L, Sancar A, Jiang Y. Genome-wide circadian rhythm detection methods: systematic evaluations and practical guidelines. Brief Bioinform. 2021. https://doi.org/10.1093/bib/bbaa135.
Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nat Methods. 2012;9:357–9.
Zhang Y, Liu T, Meyer CA, et al. Model-based analysis of ChIP-Seq (MACS). Genome Biol. 2008;9:R137.
Quinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics. 2010;26:841–2.
Ori APS, Bot MHM, Molenhuis RT, Olde Loohuis LM, Ophoff RA. A Longitudinal Model of Human Neuronal Differentiation for Functional Investigation of Schizophrenia Polygenic Risk. Biol Psychiatry. 2019;85:544–53.
Gazal S, Finucane HK, Furlotte NA, et al. Linkage disequilibrium-dependent architecture of human complex traits shows action of negative selection. Nat Genet. 2017;49:1421–7.

No competing interests reported.

Download PDF

Journal Publication

published 01 Oct, 2022

Read the published version in European Neuropsychopharmacology →

Version 1

posted

You are reading this latest preprint version

Fibroblasts as an in vitro model of circadian genetic and genomic studies

Status:

Journal Publication

Version 1

Abstract

Figures

BACKGROUND

RESULTS

Temporal RNA-seq captures genes with distinct longitudinal expression patterns

Stratified Linkage Disequilibrium Score Regression (sLDSC) analysis

DISCUSSION

MATERIALS AND METHODS

CONCLUSIONS

Declarations

References

Additional Declarations

Supplementary Files

Status:

Journal Publication

Version 1