Evidence for ~12-h ultradian gene programs in humans

Human participants and study protocol

The study protocol was approved by the University of Pittsburgh Institutional Review Board (Study 20020034; approval date: 6/4/2020) and written consent was obtained from all study participants. We studied 3 healthy male participants, who were recruited through online advertisements. Inclusion criteria consisted of individuals 18–35 years of age with a self-reported regular nighttime sleep schedule and a body mass index between 18.5–24.9 kg/m². Volunteers were excluded if they admitted to nighttime shift work or other regular nighttime sleep-disrupting activities, if they had any chronic medical conditions, took any medications or recreational drugs, or used tobacco products. Potential volunteers presented for a screening visit, inclusive of measurement of body weight, height, BMI, and laboratory studies, including a comprehensive metabolic panel (electrolytes, kidney function and liver function tests), complete blood count, 25-OH vitamin D, and thyroid stimulating hormone level, to screen for potential subclinical chronic diseases. We excluded participants with low hemoglobin/hematocrit, abnormal thyroid function and individuals with 25-OH vitamin D < 20 ng/mL.

Qualifying study participants maintained a food diary, which was used to estimate their daily caloric intake and subsequently presented to the University of Pittsburgh Medical Center (UPMC) Clinical Translational Research Center (CTRC) for a 3-day inpatient visit. On the morning of admission, participants selected items from a food menu designed to match their standard daily caloric intake with a uniform macronutrient composition of 55% carbohydrates, 25% fat, 20% protein per day. No interventions were performed during the first 24-h period of acclimatization to the hospital. Each study participant was housed in the same room in the CTRC, which contained a window to the outside. The participants were not entrained to any specific environmental cues during the study. They were encouraged to maintain their usual sleep-wake cycles; however, they maintained control over the lighting in the room. They were free to ambulate in the hallway of the CTRC; however, they did not leave the unit during the 3-day stay. The three meals were delivered at roughly the same time each day based on the standard CTRC schedule (8 A.M., 12 P.M., 6 P.M.); however, there were no rules dictating the duration of their meals. These measures were implemented to mimic free-living conditions. On the morning of the second day, an intravenous (IV) line was placed to commence blood collections at 8 A.M. and then every 2 h for 48 h (total 24 samples). Nighttime blood collection was performed through a long IV line from outside the room so that the participant would not be exposed to light or disturbed during blood collection. Blood was immediately processed in the Center for Human Integrative Physiology, two floors above the UPMC CTRC by a rotating study team, all of whom were trained in the processing procedures for this study. Blood was centrifuged (1900 RCF x 10 min) and the buffy coats were collected and immediately snapped frozen in liquid nitrogen for storage at −80C.

RNA-Seq and data analysis

RNA was isolated from peripheral blood buffy coat samples on the automated Chemagic 360 (Perkin Elmer) instrument according to the manufacturer’s instructions. Extracted RNA was quantitated by Qubit™ RNA BR Assay Kit (Thermo Fisher Scientific) followed by RNA quality check using Fragment Analyzer (Agilent). For each sample, RNA libraries were prepared from 100 ng RNA using the KAPA RNA HyperPrep Kit with RiboErase (Kapa Biosystems) according to manufacturer’s protocol, followed by quality check using Fragment Analyzer (Agilent) and quantification by qPCR (Kapa qPCR quantification kit, Kapa biosystem) on the LightCycler 480 (Roche). The libraries were normalized and pooled, and then sequenced using NovaSeq6000 platform (Illumina) to an average of 40 M 101 bp paired end reads per sample. Low-quality reads and adapter sequences were trimmed from the raw sequencing data with Trimmomatic⁶². The remaining reads were then aligned to human reference genome hg38 with STAR aligner⁶³. Gene counts were quantified with the STAR-quantMode GeneCounts function. Fragments per kilobase of transcript per million mapped fragments (FPKM) were quantified with Cufflinks⁶⁴.

Identification of the oscillating transcriptome

Averaged FPKM values at each time were used for cycling transcripts identification. Lowly expressed transcripts were removed by calculating the background expression in each participant using the average expression of a panel of 62 genes known not to be expressed in peripheral blood cells (Table S2). Temporal transcriptomes were subjected to linear detrend prior to identification of oscillations by either the eigenvalue/pencil or RAIN methods. For the eigenvalue/pencil method^11,13, a maximum of four oscillations were identified for each gene. Criterion for ~24-h genes were: period between 20 h to 25 h for first and second participants and 24 h to 30 h for the third participant, decay rate between 0.8 and 1.2; for ~12-h genes: period between 9.6 h to 13.6 h for the second and third participants and 10 h to 14 h for the first participant, decay rate between 0.8 and 1.2; for ~8 h genes: period between 6 h to 8 h for the first participant and 7 h to 9 h for the second participant, decay rate between 0.8 and 1.2; for ~16 h genes; period between 14 h to 18 h for the third participant. The relative amplitude was calculated by dividing the amplitude by the mean expression value for each gene. To determine FDR, we used a permutation-based method that randomly shuffles the time label of gene expression data and subjected each permutation dataset to the eigenvalue/pencil method applied with the same criterion⁶⁵. These permutation tests were run 5000 times, and FDR was estimated by taking the ratio between the mean number of rhythmic profiles identified in the permutated samples (false positives) and the number of rhythmic profiles identified in the original data. Analyses were performed in MatlabR2017A. RAIN analysis was performed as previously described in Bioconductor (3.4) (http://www.bioconductor.org/packages/release/bioc/html/rain.html) with either 48-h continuous data or 24-h data with biological duplicates as input¹². Genes exhibiting a period range between 10-h and 14-h and a period range between 22-h and 26-h are considered as ~12-h and ~24-h genes, respectively. FDR was calculated using the Benjamini-Hochberg procedure. Heat maps were generated with Gene Cluster 3.0 and TreeView 3.0 alpha 3.0 using Z score normalized values.

For meta-analysis, we used Fisher’s method, which combines extreme value probabilities from each test, commonly known as “p values”, into one test statistic (X2) using the formula

where pi is the p value for the ith hypothesis test. For meta-analysis of ~12-h (10–14-h) and ~24-h genes (22–26-h), the smallest p value for each gene within these period ranges was used for each individual. For example, XBP1 was found to oscillate with periods of 12-h (p = 0.068), 12-h (p = 0.014) and 10-h (p = 0.049) in participant 1, 2, and 3, respectively, and the meta p values for being a common ~12-h gene (cycling with a period between 10-h and 14-h) were calculated by “merging” the three respective p values (meta p = 0.0029). The same procedure was also performed for genes cycling with a period between 22-h and 26-h. This meta-analysis is feasible because the statistical test for each gene in each individual shares the same null hypothesis: absence of rhythms.

For RAIN analysis on temporal IR events and splicing gene expression, raw data was subjected to polynomial detrend (n = 2) before RAIN analysis.

Defining oscillating genes

The eigenvalue method can detect multiple superimposed oscillations. Therefore, we defined a ~ 24-h gene as one that exhibited a ~ 24-h rhythm, regardless of its amplitude relative to other superimposed oscillations. Similar criteria were applied to other oscillations. As such, a gene can meet criteria for both a ~ 24-h and ~12-h gene. By comparison, we define a dominant ~24-h gene as one in which the superimposed ~24-h rhythm has the largest amplitude among all oscillations. With this definition, dominant ~24-h and dominant ~12-h genes are mutually exclusive.

Intron retention detection

Intron retention events were detected by tool iREAD⁴⁵. Intron retention events are selected either with default settings T > = 20, J > = 1, FPKM > = 2 and NE score > =0.9 or more stringent settings where T > = 20, J > = 1, FPKM > = 3 and NE score > =0.9.

Gene ontology analysis

DAVID (Version 2021)⁶⁶ (https://david.ncifcrf.gov) was used to perform Gene Ontology analyses. Briefly, gene names were first converted to DAVID-recognizable IDs using Gene Accession Conversion Tool. The updated gene list was then subject to GO analysis using all Homo Sapiens as background and with Functional Annotation Chart function. GO_BP_DIRECT, KEGG_PATHWAY or UP_KW_BIOLOGICAL_PROCESS were used as GO categories. Only GO terms with a p value less than 0.05 were included for further analysis.

Motif analysis

Motif analysis was performed with the SeqPos motif tool (version 0.590) embedded in Galaxy Cistrome using all motifs within the homo sapiens reference genome hg19 as background. LISA analysis was performed using web tool (http://lisa.cistrome.org/).