Collecting mouse livers for transcriptome analysis of daily rhythms

Summary Molecular daily rhythms can be captured by precisely timed tissue harvests from groups of animals. This protocol will allow the investigator to identify transcriptional rhythms in the mouse liver while also providing a template for similar analyses in other whole metabolic organs. We describe steps for mouse entrainment, liver dissection, and rhythmicity analysis from total RNA sequencing data. The resulting rhythmic transcriptome will provide the user with a starting point for defining specific biological processes that undergo daily rhythms. For complete details on the use and execution of this protocol, please refer to Koronowski et al. (2019). A similar protocol for interfollicular epidermal cells is demonstrated in Welz et al. (2019).

darkness. This protocol can be adapted to identify purely circadian rhythms by maintaining the mice in constant darkness for two days prior to tissue harvest (Hughes et al., 2017).
Optional: Feeding behaviour has a circadian rhythm in mice;~75% of food is consumed during night-time, in the active phase (Atger et al., 2017). This feeding rhythm is a zeitgeber (time giver) or synchronizing cue for many metabolic tissues (Weger et al., 2021). To interpret results of this protocol, it is important to determine if your mutant mouse or experimental condition exhibits markedly altered feeding behaviour (as is the case in clock-disrupted animals) (Lamia et al., 2008). To eliminate food intake as a variable, restrict food access exclusively to the night-time for four days prior to tissue dissection (Atger et al., 2015). troubleshooting 1 2. Allow mice to entrain (synchronize) to the light schedules for two weeks. Change cage bedding on day seven if necessary.
CRITICAL: Avoid entering rooms during the dark period if it exposes mice to light even briefly. Ensure any light-emitting devices in the entrainment rooms are covered. Even very low intensity light (~5 lux) can perturb certain metabolic rhythms (Borniger et al., 2014;Cissé et al., 2017). If necessary, dim red headlamps can be used to check on mice during the dark period. Night vision goggles equipped with an infrared beam are ideal, finances permitting (Jud et al., 2005).
CRITICAL: Avoid creating potentially stressful events for the mice, such as loud noises, strong odors (colognes and perfumes), or bedding changes close to the date of dissection (Sorge et al., 2014). Circadian rhythms are tied to stress physiology. For example, glucocorticoids are potent zeitgebers or synchronizing cues for circadian clocks (Dickmeis, 2009;Oster, 2020).

KEY RESOURCE TABLE
Note that RNA extraction and sequencing reagents/protocols can be adapted as preferred. REAGENT  1. Take mice cages one at a time to a dedicated dissection/collection room ( Figure 2). a. Euthanize mouse by CO 2 asphyxiation, cervical dislocation, or another approved method. For dark-phase time-points, euthanize mice in darkness.
Timing: 60 s b. Decapitate mouse as a secondary confirmation. This is an important step when collecting dark-phase time-points.
Timing: 10 s c. With the mouse on its back, take small scissors and remove the layer of skin up the midline (Figure 2A).
Timing: 30 s d. Make a second cut of the muscle to expose the abdominal cavity and organs ( Figure 2B).
Timing: 30 s e. Cut down both sides, at the ribcage, to gain access to the entire liver ( Figure 2C). f. Identify the different lobes of the liver. While gently holding the left lateral lobe with forceps, cut a small piece of liver (~0.5 cm, enough for RNA extraction) ( Figure 2C).
Timing: 1 min g. Place liver piece in a 2 mL tube and immediately snap-freeze in liquid nitrogen ( Figure 2E).
Timing: 30 s h. Remove the remainder of the liver, taking care to avoid any non-liver tissue and connections.
For example, while taking care to not spill its contents, remove the gallbladder and its connections with the liver. Wrap the liver in aluminum foil ( Figure 2E). Immediately snap-freeze foil in liquid nitrogen. Pause point: Livers can be stored long-term at -80 C.
CRITICAL: When collecting dark-phase time-points, do not expose mice to light at any point. While transporting mice to the dissection room, they can be covered with heavy black trash bags or a blackout curtain. Do not turn on dissection room lights until mice have been decapitated. If collecting brain, remove the eyes before turning on the lights. CRITICAL: Ensure collection is complete within G30 min of the time-point. Consider the number of mice, data gathering (e.g., blood glucose, body weight), and tissue dissection time. For instance, if there are 10 mice that each take 6-min to process, start the harvest at 30 min prior to the time-point in order to finish at 30 min post-time-point. If possible, share the work between several investigators to reduce the amount of time spent per mouse.
Optional: In addition to liver, collect other tissues of interest. Be sure to factor in extra dissection time. Consider beforehand the downstream assays for those tissues and apply appropriate collection methods. Timing: 2 days 2. Extract RNA from the piece of liver in the 2 mL tube using an RNeasy plus mini kit (Qiagen), in accordance with the manufacture rs instructions (https://www.qiagen.com/us/shop/automatedsolutions/sequencers/rneasy-plus-mini-kit/) CRITICAL: Before proceeding to sequencing, validity of the collection should be checked by qPCR for clock genes (e.g., Bmal1, Dbp, Cry1, Per2, Rorc, Reverba). High-amplitude rhythms with distinct phases should be apparent in wild-type mice ( Figure 3; Table S1).
troubleshooting 2 3. Prepare libraries for sequencing, in accordance with the manufacture rs instructions of a total RNA library preparation kit (https://www.illumina.com/products/by-type/sequencing-kits/library-prepkits/truseq-rna-v2.html), and subsequently sequence the libraries. Sufficient depth is~20 million reads/sample for single-end sequencing.
Note: The type of sequencing, and library choice, should be directed by the hypothesis the researcher intends to test. Both paired-end sequencing and libraries enriched for specific RNA types (e.g., polyA-enriched) are compatible with this protocol.
4. Process the resulting sequencing data through a standard RNA-seq pipeline to quantify gene expression. FPKM, RPKM, and/or TPM values are compatible with downstream rhythmicity analysis.
CRITICAL: High-quality, correctly processed RNA-seq data is indispensable for identification of rhythmic transcripts using this protocol. Users are referred to bioinformaticscore-shared-training.github.io/RNAseq-R/ and bioconductor.org/packages/devel/workflows/

OPEN ACCESS
vignettes/rnaseqGene/inst/doc/rnaseqGene.html as references for standard good practice in quality control and processing of RNA-seq data.
Note: Apply Principal Component Analysis, Hierarchical Clustering or more specialised tests (see: Chen et al., 2020;George et al., 2015) to the RNA-seq data, to identify outliers that could potentially compromise rhythmicity analyses.

EXPECTED OUTCOMES
Success in this protocol will yield high-amplitude rhythms of clock gene expression with distinct phases (Figure 3), as well as rhythms in hundreds to thousands of other transcripts . As mentioned above, this protocol is designed to identify daily rhythms (see: 'limitations'). The circadian clock drives many of these rhythms, however this protocol cannot distinguish precisely the drivers of each transcript. Some transcript rhythms are not resulting from direct clockoutput genes, but may be driven by other factors such as the light-dark cycle (Li et al., 2020;Thaben and Westermark, 2016) or feeding-fasting cycle (Greenwell et al., 2019;Vollmers et al., 2009;Weger et al., 2021). Moreover, recent studies suggest that a number of transcripts can oscillate independently of the canonical core clock machinery in ex vivo and in vitro systems (Ray et al., 2020). Thus, indirect and possibly core-clock-independent daily transcripts may also be identified by this protocol. Additional experiments need to be performed to determine if a particular rhythm is circadian (constant darkness) or a direct clock output gene (disrupted in tissue-specific clock mutant mice).

QUANTIFICATION AND STATISTICAL ANALYSIS Daily rhythmicity analysis
Timing: 1 day Here we provide steps to identify transcripts with daily rhythmicity in whole transcriptome data, using the R package JTK_CYCLE (Hughes et al., 2010). Although basic knowledge of the R programming language is required, running JTK_CYCLE should not pose an obstacle even to first-time users. Many excellent R tutorials for beginners are available for free, including through The R Project: https://cran.r-project.org/manuals.html. Often, researchers will want to perform differential rhythmicity analyses to identify how the rhythmic transcriptome changes after genetic, environmental or pharmacological interventions. A discussion of how this can be achieved can be found in the 'limitations' section of the protocol.
1. Install the R programming environment. a. R can be downloaded at: https://cran.r-project.org/ 2. Download the scripts required to run JTK_CYCLE and move them to an appropriately named directory. a. The folder containing these scripts and a detailed user guide can be found at: https:// openwetware.org/wiki/File:JTKversion3.zip Recommended: The manuscript and user guide written for JTK_CYCLE by its creators describe in detail the underlying mathematics and give thorough technical advice, respectively. We highly recommended reading both of these informative publications before proceeding.
3. Using the previously-generated circadian RNA-seq data, prepare the annotation and data matrix files that will serve as the input to JTK_CYCLE (Figures 4A and 4B). a. The annotation file should comprise a first column listing the IDs of all genes quantified by RNA-seq, followed by columns containing other gene annotations (HUGO symbol, Ensembl ID, etc.) that the user prefers ( Figure 4A). CRITICAL: Ensure that the order of gene IDs is identical in the first columns of the data matrix and the annotation files.
4. In the same folder containing the scripts required to run JTK_CYCLE, save both the data matrix and the annotation files as tab-delimited text files (.txt). 5. In a text editor, open the ''Run_JTK_CYCLE (Example1).R'' script file.
a. This file can be found inside the folder downloaded from the JTK_CYCLE webpage.
6. Modify the following terms in the script to reflect the dataset being provided to the algorithm: a. On line four, alter the project name to an appropriate description (default: "Example1"). b. On lines seven and eight, change the file names to those of the annotation and the data matrix files, respectively (default: "Example1_annot.txt" and "Example1_data.txt"). (C) A "JTK_CYCLE Output File" generated from the example input files in (A) and (B), and sorted by the statistical significance of circadian gene expression. The first three columns list the different gene annotations present in the "annotation file", and the following five columns give the circadian parameters of each gene as calculated by JTK_CYCLE. In the final columns, the RPKM values used by the algorithm are repeated. The first 10 rows and 10 columns of data are shown, along with the column headers. In this example, JTK_CYCLE identified 2035 rhythmic transcripts (ADJ.P <0.01) in WT liver, which represented 11.25% of the total measured transcripts. BH.Q, Benjamini-Hochberg q-value; ADJ.P, adjusted p-value; PER, period; LAG, phase; AMP, amplitude.

OPEN ACCESS
c. For this example, an experiment with six time-points collected over 24 hr and four replicates per time-point, the following lines of the script must also be altered to reflect the dataset structure: i. Line 12, "jtkdist(6,4)" indicates the number of time-points (six) and the number of replicates provided for each time-point (four).
Note: If group sizes are unequal due to excluded or missing replicates, an additional line of code must be added to specify the number of replicates at each time-point. For example, if only three replicates are available for the first time-point, the line "group.size <-c(3,4,4,4,4,4)" would be added, and the term "jtkdist" changed to "jtkdist(length(group.sizes),group.sizes)". ii. Line 14, "periods <-6:6" specifies the range of circadian periods that are of interest. In this instance, only genes with a 24-hr period are of interest; thus, "6:6" is the input (six timepoints per cycle). iii. Line 15, "jtk.init(periods,4)" indicates the separation in hours (four) between time-points. 7. Change the working directory of R to the directory containing the data matrix file, annotation file, and scripts required to run JTK_CYCLE. 8. Paste the modified code into the R terminal and run the script.
a. JTK_CYCLE will output a tab-delimited text file containing the results into the same directory ( Figure 4C). This file provides the calculated Benjamini-Hochberg q-value (BH.Q), adjusted pvalue (ADJ.P), period (PER), phase (LAG), and amplitude (AMP) for all genes analyzed.
CRITICAL: The JTK_CYCLE script must be modified correctly to accurately reflect the structure of the RNA-seq data generated.
Note: JTK_CYCLE will typically take 10-15 min to run using a standard desktop computer or laptop.
9. To obtain a final list of robustly-rhythmic genes, the output file must be filtered to retain only those genes with clear daily rhythms in expression. Typically, an adjusted p-value of 0.01 provides a sufficiently stringent cut-off that retains only genes with clear, biologically-relevant rhythmic expression. In addition, selecting genes by expression amplitude can complement a p-value filter, adding an additional level of stringency to the rhythmic transcriptome identified.
Note: The high stringency of JTK often leads to large numbers of false negatives if the computed q-value is used for filtering (Hutchison et al., 2015). In our experience, the adjusted p-value calculated by JTK serves as a more appropriate filter for identifying genes with reproducible rhythmic expression.
CRITICAL: The p-value filter selected must be directed by the level of noise intrinsic to the dataset obtained, which will vary on a case-by-case basis. Visual inspection of core clock genes, such as Bmal1, Per2 and Cry2, can provide a good starting point to assess the technical and biological noise present in a circadian dataset (Figure 3).
10. The final set of circadian genes can now be used for further downstream bioinformatics analyses, such as gene ontology and pathway analysis, or to direct future in vivo and in vitro experiments.

Note:
For examples of phase plots, amplitude plots, and expression heat maps typically used to visualize circadian gene sets, please refer to Welz et al., 2019Welz et al., 2019).

LIMITATIONS
Here we describe how to collect liver tissue over one circadian cycle, with a large number of replicates per time-point and a relatively long sampling interval (4 hr). The experimental design used to characterize the ll OPEN ACCESS STAR Protocols 2, 100539, June 18, 2021 daily rhythmic transcriptome varies extensively across published circadian studies. As such, the design employed to answer each new research question should be determined by balancing cost versus time considerations, and considering the particular focus of the study. An alternative, and arguably more stringent, strategy to that employed in this protocol would be to collect samples over two consecutive circadian cycles (48 hr) and to reduce the sampling interval (e.g., collect every 3 hr). This would make the data more robust to outliers and reduce false-negatives in rhythm detection (Hughes et al., 2017). However, due to cost and time limitations, the number of biological replicates is often sacrificed in such studies. This often makes it more difficult to calculate circadian parameters that require large numbers of biological replicates to provide the appropriate statistical power and precision. A detailed set of established guidelines for circadian 'omics' experiments is provided by Hughes et al., 2017, and it is highly recommended that the user study these guidelines in some detail before deciding on their final experimental design (Hughes et al., 2017).
The protocol described here requires liver sampling from different mice at multiple time-points. Currently, it is not possible to perform repeated measurements of the liver transcriptome in the same mouse over the circadian cycle. Thus, this experimental strategy has its limitations when obtaining samples and interpreting data from mice that may lack phase-alignment, or if the amplitude of the circadian rhythms in the individual experimental animals varies (e.g., after prolonged free running conditions in complete darkness). Under such experimental conditions, alternative methods might have to be applied to determine circadian rhythmicity. For example, clock gene rhythms can be visualized by sequential ear snip biopsies from the same mouse over a circadian cycle .
In this protocol we have described a procedure for identifying the daily rhythmic transcriptome under physiological light-dark conditions. To formally demonstrate that the rhythmic transcription identified is indeed driven by the circadian clock, or oscillates in a circadian manner independently of the core clock, it must meet a number of established criteria (Dunlap, 1999). First, the free-running period of the oscillation must be temperature compensated i.e., remains~24 h across a broad range of physiological temperatures. Second, the oscillation must be 'entrainable' to external cues i.e., the phase of the gene shifts in response to changing environmental conditions, such as light cycles. Finally, a circadian oscillator must be self-sustained i.e., continuing to oscillate under constant environmental conditions, such as continuous darkness. A full description of the established methods used to test these conditions is outside the scope of this protocol, but users are referred to Eckel-Mahan and Sassone-Corsi as a reference for such approaches (Eckel-Mahan and Sassone-Corsi, 2015).
A wealth of algorithms has been developed to identify circadian features of "omics" data, each of which has specific intentions, strengths, and weaknesses. Although JTK_CYCLE is one of the most established and frequently cited examples for rhythmicity detection, we recommend validating results using different algorithms (Guan et al., 2020;Koronowski et al., 2019;Welz et al., 2019). Specifically, users are directed to BIO_CYCLE and MetaCycle as complementary approaches for rhythmicity detection (Agostinelli et al., 2016;Wu et al., 2016). BIO_CYCLE applies machine learning techniques, operates as a web tool, and provides a rapid and effective alternative to JTK. MetaCycle, in contrast, is an R package that enables the user to apply three popular rhythmicity detection algorithms to the same data set (ARSER, JTK_Cycle and Lomb-Scargle), allowing for independent confirmation of rhythmicity and avoidance of false-negatives. Moreover, when comparing the rhythmic transcriptomes of multiple groups (genotypes or conditions) users should consider employing algorithms specifically designed for differential rhythmicity analysis. The R package DryR (Differential Rhythmicity Analysis in R) enables a combined rhythmicity detection and differential rhythmicity analysis, allowing the user to detect changes in rhythmicity between multiple groups (Weger et al., 2021).

Potential solution
Holding/breeding rooms can be used for entrainment. However, disruptions from staff or other researchers may compromise circadian behavior and thus will present a major limitation to the experiment. Alternatively, a range of cage and cabinet solutions have been developed to overcome such space limitations (cage -www.tecniplast.it/uk/product/leddy-lighting-control-system-with-in-gm500red-cage.html, cabinet -www.actimetrics.com/products/circadian-cabinets/). These cages and cabinets enable the investigator to subject mice to alternative light schedules without the need for specific rooms, whilst also minimizing unwanted disruption to the mouse's normal daily rhythms.

Problem 2
Clock gene expression is non-rhythmic or irregular (step 2 in ''step-by-step methods details'').

Potential solution
Check that entrainment parameters are within the desired range. Measure humidity and temperature, and confirm that the room lights are turned on and off at the correct times using light meters able to log patterns over multiple days. Monitor entrainment rooms for disturbances during the day (for instance, other laboratory members spending time in the rooms, or entering and exiting the rooms frequently). Notify colleagues of entrainment weeks; hang signs on doors to indicate that there is an ongoing, highly sensitive experiment. During entrainment, regularly ask the facilities management about any major disruptions in power that might have occurred.

Problem 3
If not enough age-and sex-matched mice can be generated at once (e.g., when using mice with a genotype of low mendelian ratio), sample collection over two days, as described in this protocol, might be impossible (step 1 in ''before you begin'').

Potential solution
When breeding mutant mice, it may not be possible to collect all samples over two days as described in this protocol. As long as experimental conditions can be kept constant, tissues can be collected over a period of months, and in our hands, this is suitable to detect rhythms by RNA sequencing. This approach may not be appropriate for other applications, such as measurement of certain unstable metabolites. When collecting over months, randomize genotypes and time-points to limit batch effects, and collect livers continually as mice become available (e.g., when reaching the right age).

Problem 4
JTK_Cycle calculates few or no genes to have daily rhythms of transcription (step 9 in ''quantification and statistical analysis'').

Potential solution
Typically, an investigator should expect to detect 1000+ genes oscillating in a daily manner when this protocol is applied to the livers of wild-type mice. Nonetheless, this number will vary substantially depending on the whether an animal has been genetically modified or subjected to some environmental/pharmacological perturbation (see: Welz et al., 2019 andKoronowski et al., 2019 for examples of the scale of this variation upon genetic modification). Aside from true variation in the size of the rhythmic transcriptome, the degree of biological and technical variation present in the RNAseq data gathered will determine the sensitivity with which daily rhythms can be identified by JTK_Cycle. Thus, it is essential that investigators: 1) Use a number of biological replicates that reflects the anticipated biological variability between mice, 2) Perform stringent quality control of the RNA-seq data input to JTK_Cycle, 3) Adjust the stringency of the filters used to classify daily rhythms according to the variability present in the data.

OPEN ACCESS
STAR Protocols 2, 100539, June 18, 2021 Problem 5 JTK_Cycle fails to execute correctly, and the expected output files are not generated (steps 8 and 9 in ''quantification and statistical analysis'').

Potential solution
Ensure that the input files have been correctly formatted in line with the examples provided in Figure 4. In addition, check that the modifications required to the JTK_Cycle script have been made correctly, without syntax errors, and that the modifications correctly reflect the structure of the RNA-seq data. If the problem persists, users should consult the user manual of JTK_Cycle or contact its developers directly.

RESOURCE AVAILABILITY
Lead contact Requests for further information and/or resources and reagents should be directed to Kevin Koronowski (Kkoronow@uci.edu).

Materials availability
This study did not generate new unique reagents.

Data and code availability
Datasets associated with this protocol are available as supplemental tables in Koronowski et al., 2019 and the accession number for the sequencing data is GEO: GSE117134 . Please also see 'key resources table' for analysis algorithm information.