The importance of circadian regulation in both chemotherapy  and radiation therapy is increasingly well documented . Conducting cancer therapy at optimal times of day based on circadian variations of sensitivity to cancer treatment (chronotherapy) has great potential for improving therapeutic efficacy and/or decreasing side effects [1, 2]. Oral mucositis (OM) is a common debilitating complication of cancer chemo- and radiotherapy occurring in nearly all patients with head-and-neck cancer and causes significant morbidity and reduced therapeutic effectiveness [3, 4]. With limited prevention and treatment options, there is an urgent need to develop novel strategies for preventing or mitigating OM to improve patients’ outcomes and quality of life . Treatment-induced OM is typically seen as an “outside-in” process, in which chemo- or radiotherapy nonspecifically targets the rapidly proliferating cells of the basal epithelium, causing the loss of ability to self-renew, inflammation and ulceration . We and others have shown that the severity of OM may depend on time of radiotherapy [2, 6, 7]: patients treated in the morning have less severe OM. These findings are biologically plausible, as many key regulators of fundamental biological processes that influence tissue response to cancer therapy, such as cell cycle progression [8, 9, 10, 11] and DNA damage response , are under circadian control . Therefore, choosing treatment times associated with less radio-sensitivity in non-cancerous oral mucosa could be one option for reducing severity of radiation-induced OM.
The application of chronotherapy is complicated by inter-individual variability of circadian phase in humans and differences between central (i.e., in the suprachiasmatic nucleus) and peripheral (eg., oral mucosal tissue) circadian phase. The Dim Light Melatonin Onset (DLMO) is currently considered to be the most accurate marker of the circadian phase in humans . However, there are two problems with using it as a marker of circadian phase for chronotherapy: (i) determining DLMO requires multiple consecutive samples of blood or saliva collected across an 8- to 24-hour interval under controlled conditions and (ii) DLMO represents the central clock time (“phase”) , which differs from the circadian clock time in peripheral tissues  that are vulnerable to the damage of cancer therapy (both chemotherapy and radiation).
Critical determinant of sensitivity of normal cells to radiation therapy is the cell cycle phase: cells in the G1 and S phases are less radio-sensitive than those in the G2/M phases [16, 17]. Circadian rhythms in mitotic index and DNA synthesis, the major events in the M and S phase of a cell cycle, have been documented in epithelium of mouse alimentary tract (intestine, tongue, esophagus and stomach) [18, 19], as well as in human rectal and oral mucosa [20, 21]; these findings provide a foundation for using markers of cell cycle progression as potential biomarkers for development of personalized chronotherapeutic schedules for cancer treatment. Previous study analyzing repeated oral mucosal samples from six healthy male volunteers has reported diurnal variation of five cell-cycle-associated proteins (Cyclin E, A, B1, Ki-67 and p53) . However, the sample collection method (invasive punch biopsy) limited the potential for future clinical application.
In this pilot study, we assessed diurnal variation in mRNA expression of genes relevant to treatment response using oral mucosa samples collected by a non-invasive brush biopsy, in healthy volunteers . Factors that may impact inter-individual difference of circadian timing were also examined. This study was designed for examining diurnal rhythm (i.e., endogenous circadian plus potential driven effects from behavior cycles such as the sleep-wake and feeding cycle), rather than endogenous circadian rhythm (purely driven by endogenous pacemaker); the latter requires a burdensome constant routine protocol to eliminate all periodic changes in behavior and environment . By allowing the participants to follow their usual schedule of sleep and meal will generate the data that reflex participants’ real-life experience, and therefore has better potential for translation to clinical setting.
This work is a fundamental first step in developing biomarkers for determining an optimal timing of cancer treatment for an individual patient to reduce oral mucositis by establishing whether the expression of these genes is rhythmic across the 24-h day.
We recruited 11 healthy volunteers, 3 men aged 23–50 and 8 women aged 24–60, with no history of cancer diagnosis, through the University at Buffalo (UB) Clinical and Translational Science Institute (CTRC) Community Recruitment Liaison, and an email circulated to UB students. Exclusion criteria for the study included: pregnant females or females using oral contraceptives or hormonal therapy; presence of open blisters, sores or lesions in the mouth that would make brushing of the buccal mucosa difficult or painful; current self-reported sleep disorders (e.g., sleep apnea, delayed sleep phase disorder, advanced sleep phase disorder); psychiatric conditions (i.e., depression or anxiety); color blindness (as determined by the Ishihara test) ; current smoking (e.g., every day or some days); drug use (e.g., marijuana, opioids, cocaine, amphetamine); inability to stop drinking coffee or alcohol for 24 hours before and during sample collection; use of any medications that may influence sleep or circadian rhythm within three weeks before and during the sample collection; working the night shift within at least two months prior to the study; and traveling across more than one time zone in the month preceding the study. Participants provided written informed consent prior to their participation. The study was approved by the Roswell Park Comprehensive Cancer Center Institutional Review Board. This study is not a clinical trial, therefore is not required to be registered on Clinicaltrials.gov.
Each participant completed the Pittsburgh Sleep Quality Index questionnaire  and the Morning-Evening Questionnaire  at consent and was asked to wear a wrist actigraphy monitor (Actiwatch Spectrum; Philips Respironics, Oregon, USA) on the non-dominant arm. This actigraphy monitor, together with a sleep log, were used to monitor the participant’s sleep-wake cycle for one week prior to sample collection . The time interval between actigraphy sleep data collection and laboratory sample collection ranged between 2 days to 32 days (<7 days for 5 participants, 11 days for one, 14 days for three, 21 days for one, and 32 days for one participant).
Volunteers were asked to arrive at the CTRC between 9:00 and 9:30 on the scheduled date. They stayed in a private CTRC room for sample collection for approximately 21 hours until about 7:00 next day; each private room had an examining table and chair, a reclining chair, a hospital bed for overnight stays, and windows with light filtering curtains. A detailed sample collection protocol is summarized in Figure 1. Briefly, six samples of full layers of cheek buccal cells were collected by a trained staff member every 4 hours, from about 10:00 to 6:00 the next day using the OralCDx® brush (CDx Diagnostics, NY), a non-invasive method for collecting cells from all epithelial layers, typically used in dental offices for early detection of pre-cancer . In parallel, saliva samples were collected every hour, starting from about 6 hours before the participant’s usual bedtime until 2 hours after their usual bedtime. Saliva was collected before cheek cell samples when collection times of saliva and cheek overlap each other. Throughout this time interval, samples were collected under dim light conditions (<5 lux). Saliva samples were collected using Salivette cotton (Sarstedt, Newton, NC) that participants were instructed to put into their mouths and chew for about 1 min until saturated, before spitting it back into the Salivette. No saliva sample was contaminated by blood, as determined by checking the color of saliva. Food was not allowed within 15 minutes before saliva sample collection. After each meal or snack, participants were instructed to brush their teeth using water without toothpaste. Tooth brushing, rinsing or water drinking was not allowed within 10 minutes before saliva collection. Volunteers were allowed to maintain their usual activities and take normal meals; they were encouraged to read, write, and listen to music while in the facility. During sleep time the rooms were dark; a dim light was allowed at the time of sample collection and when participants needed light in the room.
Immediately following each brush mucosal cell collection, the brush was immersed in Qiazol (cell lysis Reagent from Qiagen, USA), mixed well, and stored temporarily in a –80° C freezer at CTRC-CRC. The next day, all the completed brush mucosal samples were delivered on dry ice to Roswell Park’s DataBank and BioRepository (DBBR) Shared Resource for storage at –80° C. After sample collection was completed from all participants, they were delivered on dry ice to the Roswell Park Genomics Shared Resource for RNA extraction and RNA-seq assay. Purified RNA was prepared using the miRNeasy micro kit (Qiagen, Cat No. 217084) and sequencing libraries were prepared with the KAPA RNA HyperPrep Kit with RiboErase (HMR) (Roche, USA), from 100ng total RNA. The final RNA-seq libraries were sequenced on Illumina NovaSeq 6000 using 100 paired end sequencing.
Bioinformatics pre-processing and quality control (QC) steps were carried out by the Roswell Park Bioinformatics Shared Resources, using the established pipeline following commonly adopted practices for RNA-seq data analysis. Raw reads that passed the Illumina RTA quality filter were demultiplexed and pre-processed by using FastQC for sequencing base quality control. The reads were then mapped to the latest version of human reference genome (GRCh38) and reference transcriptome GENCODE (v25) using Bowtie (v1.0.1)  and TopHat (v2.0.13)  aligner. A second round of quality control was performed to identify problematic samples with RNA-seq library preparation using the RSeQC (v2.3.2) . Because oral mucosal samples contain bacteria, which will decrease the mRNA integrity number (mRIN), we additionally checked gene body coverage of housekeeping genes to better assess the integrity/degradation of the human RNA samples. Transcripts Per Kilobase Million (TPM) estimates for gene expression were obtained using Salmon (v0.13.1). A total of 19798 genes had detectable expression data and were qualified for further analysis. Log10 transformed TPM values with a 0.1 correction constant were used for harmonic gene testing (MetaCycle).
The RNA qualities from all samples were acceptable with over 3 million reads obtained to quantify gene expression. The mRIN ranged from 2.0 to 6.8 (Supplementary Figure 1.), acceptable for whole transcriptome (ribosomal depletion) library preparation. The mapped reads were evenly distributed onto different gene features (CDS exons, UTRs) and the coverage across whole gene body did not show obvious significant 3’ or 5’ bias for any sample (Supplementary Figure 2a.). Due to the nature of sample collection site (i.e., oral cavity), microbiome RNA contamination was inevitable, and the mapping rates reflect the amount of foreign RNA in the sample (Supplementary Figure 2b.). To better assessing the RNA integrity, Supplementary Figure 3. showed that the house keeping gene body coverages were overall uniform in all samples, which means that there was no obvious degradation of human RNA happening in any sample. So, the low RIN (<4) for some samples is primarily caused by microbiome RNA and does not reflect the sample quality. After mapping and assigning reads to human genome reference using Salmon, a mean of 10,689,349 and sd of 8,075,814 reads were quantified in total.
Saliva samples were centrifuged immediately to extract the saliva from the cotton swab and then frozen at –80° C temporarily at CTRC-CRC. The next day, all saliva samples were delivered on dry ice to Roswell-DBBR for storage at –80° C. After completing field sample collection from all participants, saliva samples were shipped on dry ice to Solidphase, Inc. (Portland, ME) for melatonin radioimmunoassay using commercially available kits (BUHLMANN, Schönenbuch, Switzerland). Based on duplicate quality control samples, the intra-assay and inter-assay CVs were 1.9% and 6.6%, respectively, for the low concentration samples (range 1.6–4.3 pg/ml); 8.0% and 0.9%, respectively, for the high concentration samples (range 16.7–31.1 pg/ml). Circadian phase was defined as the clock time at which melatonin concentration exceeded the threshold of 3 pg/mL, indicating the onset of melatonin secretion under dim light conditions [Dim Light Melatonin Onset (DLMO)] estimated by linear interpolation based on the clock time of two consecutive samples with melatonin concentration less than and greater than the threshold, respectively. Specifically, if the two known values for the adjacent times around the threshold were (t1, m1) and (t2, m2), then the time of.
To identify genes whose expression had a significant variation by time of a day, we used a cosine model implemented by the MetaCycle R package . The MetaCycle estimates circadian parameters from three independent rhythm detecting models (ARSER, JTK_CYCLE and Lomb-Scargle) and integrates them using a voting scheme. The three methods considered by MetaCycle have their own strengths and weaknesses depending on how the data was generated. In the case of the present study, the data are continuous, unevenly spaced with 6 time points per subject which better suits the features covered by the Lomb-Scargle method. Briefly, Lomb (1976), proposed a model using least-squares to fit sinusoidal curves. Scargle (1982) extended this method by defining a Lomb-Scargle periodogram and obtaining the null distribution for it . This fit is conducted on the gene expression data of six time points by each gene and participant. Resulted estimates of circadian parameters for each participant were then integrated using the MetaCycle’s meta3d function.
Meta3d can directly output the p-value of rhythmic signal, phase, amplitude, relative amplitude values (rAMP) and baseline value. Amplitude is half the difference between the estimated minimum and maximum values (23 and 25 were set for this study) Relative amplitude values are the ratio of the amplitude and baseline value, which may be used to compare the amplitude values among genes with strongly different expression levels. Group-level average phase was calculated using the circular mean from each participant’s phase estimate. For details of parameter estimation in Meta3d, refer to the link (https://cran.rproject.org/web/packages/MetaCycle/vignettes/implementation.html). We selected a priori candidate genes involved in circadian regulation, cell cycle progression (a strong predictor of radio-sensitivity), and DNA homologous recombination repair (a critical pathway for the repair of DNA double-strand breaks). We focus on these pathways because they are relevant to circadian regulation and/or the response to cancer therapy [16, 17]. The gene list included 16 well-established circadian genes with noted circadian variation in human skin and mouse tissues : ARNTL, NPAS2, CLOCK, PER1–3, CRY1–2, CIART, NR1D1, BHLHE41, DBP, TEF, HLF, RORC, NFIL3. The genes for the cell cycle progression (n = 118) and homologous recombination DNA double-strand break repair (n = 26) were based the Molecular Signature Database  using c2-CP (curated gene set – canonical pathways). According to a methodological guideline for genome-scale analysis of biological rhythms , we used different thresholds to define whether a gene mRNA expression was rhythmic. We present major findings based on the criteria of p-value ≤ 0.1 and relative amplitude (rAMP) > 0.1. False discovery rate  was then used to adjust for multiple comparisons for number of genes examined in each pathway.
To evaluate the peak time variability of these diurnally expressed genes, we calculated the range of point estimate for gene peak times, relative to the wall clock time, circadian phase and sleep time. The gene expression peak times relative to circadian phase or sleep time were calculated by subtracting the circadian (DLMO) time or the midpoint of sleep-onset and wake-up (midsleep) time from original peak times.
We further assessed the circular correlation between the peak times of identified genes and phases of other biological rhythms (the DLMO time and midsleep time) using the CircStats R package . Specifically, the “circ.cor” function was employed to estimate a circular version of Pearson’s product moment correlation. Distribution of variables were examined using goodness of fit test from the package to ensure the applicability of the correlation calculation.
To examine correlation between gene peak times and characteristics of participants (age, sleep time), we used Partial Spearman correlation tests (SAS 9.4). Each analysis adjusted for covariates of age, sex, sleep time, and DLMO, excluding the variable that was being examined.
Using gene expression data of all 66 samples from 11 participants, we examined the correlation between mRNA expression of circadian genes and WEE1, a gene involved in G2/M cell cycle checkpoint inhibition, using repeated measures correlation method accounting for intra-participant variation as implemented in R package rmcorr (v0.3.0) . Correlation values are displayed and organized using complete Euclidean clustering, and the heatmap was plotted using heatmap R package.
Demographic information about the study participants, including age, sex, race, chronotype and sleep times is summarized in Table 1. We studied 8 women and 3 men, 4 different ethnicities, ages ranging from 23 to 60 years. Based on the Morningness-Eveningness questionnaire , 9 participants belong to moderate morning or intermediate chronotype, 1 to moderate evening and 1 to definite morning type. The DLMO ranged from ~17:30 to 22:00. The PSQI-based bedtimes were between 20:30 and 1:30, wake-up times were between 5:45 and 7:45, and midsleep times (between bedtime and wake-up time) were between 1:15 and 3:15. Actigraphy-based wake-up time ranged from 5:31 to 7:57 on weekdays and from 3:19 to 8:40 on weekends. The midsleep time ranged from 1:56 to 4:07 on weekdays and from 23:28 to 6:01 on weekends.
|VARIABLES||NO.||%||MEAN (STD), MIN, MAX|
|Chronotype defined by the morning-evening questionnaire scores (Horne and Ostberg, 1976, ref 23)|
|Moderate evening: 31–41||1||9.1|
|Moderate morning: 59–69||5||45.5|
|Definite morning: 70–86||1||9.1|
|DLMO (3pg/ml) time (HH;MM)||20:13 (1:25), 17:22, 22:04|
|Bedtime-PSQI* (HH:MM)||22:42 (1:18), 20:30, 25:30|
|Waketime-PSQI* (HH:MM)||6:24 (0:36), 5:45, 7:45|
|Midsleep-PSQI* (HH:MM)||2:24 (0:42), 1:15, 3:15|
|5-day average weekday midsleep – Actigraphy** (HH:MM)||3:01 (0:48), 1:56, 4:07|
|2-day average weekend midsleep – Actigraphy** (HH:MM)||3:39 (1:46), –0:28, 6:01|
|5-day average weekday waketime – Actigraphy** (HH:MM)||6:25 (0:47), 5:31, 7:57|
|2-day average weekend waketime-Actigraphy** (HH:MM)||7:16 (1:17), 3:19, 8:40|
|Intervals between completion of actigraphy data collection and sample collection|
The number of genes that can be defined as displaying diurnal variation depended on significance thresholds. Out of 19798 genes detected, 3 genes had p-value < 0.01, 108 genes had p-value < 0.05, 650 genes had p-value < 0.1, 2800 genes had p-value < 0.2 and 5514 genes had p-value < 0.3. Among the candidate pathway genes selected a priori, the three genes with the lowest p-values, PER3, CIART, and WEE1 reached a criterion of p-value < 0.01 and relative amplitude (rAMP) > 0.1 (Table 2). After multiple comparison adjustment for number of genes tested in each pathway, PER3 and CIART remained significant (FDR < 0.1). In addition, the circadian pathway genes PER1, ARNTL, TEF, CRY2, and PER2 had p-value ≤ 0.1 and rAMP > 0.1 (Table 2). The average peak times of RNA expression were approximately 10:15 for PER3, CIART and TEF, 10:45 for PER1, followed by PER2, CRY2 and WEE1 peaking about 13:00; whereas ARNTL peaked about 19:30. Results of all candidate genes in the three pathways and the 108 genes with p-value < 0.05 are presented in Supplementary Table 1.
|GENE||AMP||RAMP||CLOCK PEAK TIME*||PEAK TIME RANGE||P-VALUE|
|POINT ESTIMATE||AVERAGE||95%CI||TIME IN CLOCK||DLMO-ADJ**||MIDSLEEP-ADJ**|
The mRNA gene expression patterns over 20 hours of data collection (i.e., 10AM – 6AM) for the most significant five genes (CIART, WEE1, PER3, PER1, and ARNTL) are presented in Figure 2. The peak times and amplitudes varied across participants (Figure 3). The range of peak time was 7.8 hours for CIART, 16.7 hours for WEE1, 8.1 hours for PER3, 8.0 hours for PER1, and 16.2 hours for ARNTL (Table 2) The range of the relative gene peak times adjusted for the DLMO and midsleep times was larger for CIART, PER3, PER1 and ARNTL and smaller for WEE1 (Table 2). The range of peak times are greatly impacted by a few outliers: removing one participant who had outlying peak times reduced the average peak times for the study sample to 8.5 hours for ARNTL (Figure 3). For PER1, eight participants had peak times between 10:30 and 14:10 and three participants had peak times at 6:30, 7:00 and 9:00.
The mRNA expression levels for PER1, PER2, PER3, CRY2, CIART, TEF, DBP and NR1D1 were positively inter-correlated with each other (p-value < 0.05) and clustered together (Heatmap in Figure 4). ARNTL was negatively correlated (r = –0.48 – –0.20) with the above-mentioned genes. WEE1 was positively correlated with PER2 (r = 0.6, p-value < 0.001), but had weaker correlation (r = 0.29–0.34, p-value < 0.05) or no correlation (p-value > 0.05) with other genes in the circadian pathway (Figure 4).
The correlations between peak times of genes with diurnal variation (Table 2), sleep time and DLMO time are shown in Figure 5. DLMO time was significantly correlated with sleep times collected by questionnaire (r = 0.68~0.79). There were also positive correlations for peak times of gene expression between CIART and PER2, CIART and PER3, PER1 and PER3, as well as ARNTL and PER2 (r = 0.61~0.66). However, peak times of genes with diurnal variation were not correlated with DLMO or sleep times, except for a negative correlation between PER1 and the DLMO time.
The peak times of PER1 and PER3 were strongly correlated with age; older people tended to have later peak times after adjusting for gender, DLMO and sleep time (Partial Spearman correlation = 0.77 and 0.69, respectively) (Table 3, Supplementary Figure 4). We also found that older age has a tendency to be associated with earlier sleep time (r = –0.38) and earlier DLMO time (r = –0.4), however the small number of samples analyzed did not allow reaching statistical significance (Supplementary Figure 4). No correlation was observed between gene peak times and sleep time (Table 3).
Using oral mucosal samples collected over 20 hours in 11 healthy volunteers by a non-invasive brush biopsy, we found that several genes involved in circadian rhythm regulation and cell cycle progression display time-of-day variations in their pattern of expression. Out of these genes, WEE1 and PER1 are of most interest, which play direct roles in cell cycle progression and apoptosis, and associate with cancer therapy response; the finding of WEE1 rhythmicity, previously reported in liver tissue of mice  and gut epithelium of rats [40, 41], has never been reported in human oral mucosa. Peak times of oscillating genes differed by individuals; such variability across individuals could at least partially be explained by age after adjusting for sex, sleep time and DLMO.
We found that the peak times for mRNA expression of PERs and CRY2 ranged from 10:00–13:30 across individuals, and ARNTL peaked on average at 19:30; there were strong positive correlations in mRNA expression of the PERs and CRY2 genes, and negative correlations with ARNTL mRNA expression. Similar peak time order and correlations of these circadian core clock genes were reported previously in human skin . These findings are consistent with the molecular mechanism of circadian rhythms, which is generated by an autoregulatory transcriptional feedback loops through positive and negative elements . In mammals, during the day, CLOCK (or NPAS2) interacts with BMAL1 to activate transcription of the PER and CRY genes, resulting in high levels of these transcripts, while during the night, the PER–CRY repressor complex is degraded, and CLOCK–BMAL1 can then activate a new cycle of transcription . Our data suggest that in addition to skin, fibroblast and hair follicles [17, 34, 44], gene expression in oral mucosal tissue in humans display strong diurnal rhythms and therefore may be considered to serve as a biospecimen for developing circadian biomarkers. This is particularly important as unlike skin and fibroblast, collection of oral mucosal tissue allows for non-invasive multiple sampling.
These results provide a biological interpretation for clinical reports associating severity of OM with delivery time of radiotherapy . In our own work using data from the electronic medical records of head and neck cancer patients, we found that the most severe cases of OM were associated with a radiation delivery time window of 11:30–15:00 , which overlaps with the peak time range of mRNA expression for several genes (Table 2) that we observed in healthy volunteers in the current study. During radiotherapy, OM is initiated by DNA damage and generation of reactive oxygen species, both of which trigger a series of interacting biological events including inflammation reactions, cell cycle deregulation, DNA damage response, and apoptosis . Of all rhythmic genes identified, WEE1 and PER1 have been of most interest for further biomarker development due to their well-reported functions in these processes. WEE1 is a tyrosine kinase that phosphorylates CDC2 and prevents cell progression through G2/M and S cell cycle checkpoint into mitosis for DNA repair in response to DNA damage . Preclinical and clinical studies have demonstrated encouraging antitumor effects of WEE1 inhibition . Therefore, we expect expression of WEE1 would be negatively associated with radiation sensitivity and radiation-induced OM. The circadian pattern of WEE1 transcription has been reported in epithelium of the rat duodenum, ileum, jejunum, and colon , and is directly activated by the CLOCK–BMAL1 transcriptional complex in mouse liver . On the other hand, ectopic expression of PER1 was found to inhibit WEE1 expression in human colon cancer cell lines . Therefore CLOCK-BMAL1 and PER1 may regulate G2/M progression via WEE1. Circadian genes were also demonstrated to affect cell cycle progression at G1/S transition. PER1 inhibits p53-induced expression of p21, which prevents G1/S transition by inhibiting cyclin-dependent kinase activity including CDK4/6 [8, 47]. These roles of PER1 on the G1/S and G2/M check points may free the cell cycle arrest from repairing DNA damage and lead to DNA damage induced apoptosis . PER1 has also been directly involved in DNA damage response by interacting with the ATM and ATR DNA damage check point kinases or modulating activity of p53 . In line with these roles in cell cycle progression and DNA damage response, the overexpression of PER1 has been shown to sensitize human cancer cells to ionizing radiation-induced apoptosis and cause significant growth reduction . Therefore, we expect PER1 expression to be positively associated with radiation sensitivity and radiation-induced oral mucositis. In the current study, 8 (out of 11) participants had PER1 expression peaking 10:30~14:10, whereas expression of WEE1, a gene expected to be negatively associated with radio sensitivity, peaked during 13:30–19:00 in 8 participants; P53 did not show obvious diurnal pattern of expression. Although the exact mechanisms might be more complicated and need further study, this peak time range of PER1 expression (10:30–14:10) largely overlaps with the delivery time window of radiotherapy (11:30–15:00) associated with the most severe OM in head and neck cancer patients, which might partly interpret the underlying mechanism of the clinical observation .
Our results confirm the hypothesis that the peak times for cycling genes display inter-individual variability, which is one of the challenges for applying chronotherapy in the clinic. The range of peak times are gene specific: some cycling gene expression levels (CIART, PER1, PER3) had smaller range (about 8 hours) than others (about 11.5 hours for PER2 and 16 hours for ARNTL and WEE1). Our results suggested such variability may be partially contributed by individual characteristics: those with older ages tended to have later peak times for PER1 and PER3. This finding, however, is not in line with the documented phase shift toward a “morning” chronotype  and a phase advance in brain tissue  with aging; similar trend was observed in our data: earlier sleep time and earlier DLMO time associated with older age. But the age effect on the circadian clock timing in peripheral tissues may differ from the age effect on central clock or brain tissue. Further study is needed to verify or negate our finding on gene peak time and age. We also found that women tended to have later peak time of circadian gene expression than men (Supplementary Table 2), but the correlation with only one gene expression (PER3) reached borderline statistical significance. Given that we only included 3 men, these results need to be interpreted cautiously and need to be confirmed in a larger study. Further studies are also needed to understand more about the role of these genes associated with radio-sensitivity. These findings may help in developing an approach to guide individual timing of chronotherapy.
We did not find correlations between gene expression peak times and chronotype, sleep times or DLMO, except for negative correlation between DLMO and PER1. The lack of correlation with circadian or sleep timing suggest that these parameters representing central clock do not directly reflect the phases of gene expression in peripheral tissues. Overall, it emphasizes the importance of assessing the circadian timing of specific peripheral tissue that is subject to side effect of cancer therapy. However, given the large confidence intervals around the peak time estimates, potential confounding factors by age, as well as the small sample size (11 participants), it is possible that associations may exist and differ from what we observed. Waiting days (2–32 days) between actigraphy sleep measurement and sample collection days also limited the capability of finding correlation between sleep timing and gene peak timing. Studies with more frequent sample collection are needed to estimate the peak times of these genes more precisely and accurately. Better design is also needed to assess actigraphy sleep time that represent sample collection days.
Our study has several limitations. As mentioned above, the small sample size of this pilot study limited the statistical power for finding genes with diurnal rhythmic expression, and for assessing correlations between gene peak times and patient characteristics. Future studies with larger sample size are needed to confirm our findings. Our analysis was based on an assumption that mRNA expression follows a sinusoidal pattern with a 24-hour period. Violation of such assumption may also influence peak time estimation. To test this hypothesis, sample collections that cover more than 24-hour time interval will be ideal but collecting more than six consecutive samples from human oral mucosa is challenging. Our sampling design (every 4-hours for 6 samples) has been used in previous studies to determine diurnal rhythmicity . Additionally, given that seven of the eight statistically significant genes belong to circadian regulation pathway and there is evidence of circadian expression of these genes in mammalian tissue, including gut epithelium [11, 32, 41, 43], the assumption that mRNA expression of these genes follow a sinusoidal pattern in oral mucosa is likely to be true. There are other factors that can introduce random variations in measurement of RNA expression in genes, including differences in cell types present within the oral cavity (including bacterial), and sensitivity of the RNA expression assays. These factors can be minimized in the future by choosing more appropriate assay protocol. In this study, a standard whole transcriptome library preparation protocol was used. As a result, bacterial RNA from the microbiome may have accounted for a substantial proportion of the total sequenced RNA, which may have reduced the quantification accuracy of low expression genes. Alternative library preparation protocols, like capture-based RNA library preparation or message RNA specific library preparation with PolyA selection, may be more suitable. A low mapping rate to human genome ratio reflects a greater microbiome population or activity. Interestingly, the microbiome population and activity in our samples suggested a circadian pattern with a lower proportion of microbiome RNA during the day and higher at night (Supplementary Figure 1., Supplementary Figure 5.), which could be further investigated. Despite the study limitations, we were able to detect diurnal variations in a number of genes in the circadian pathway and cell cycle pathway. In addition, the correlations in peak times between the core circadian genes were consistent with their biological interactions and previous reports, strengthening the validity of our findings.
In conclusion, we found diurnal variations in mRNA expression patterns for genes involved in circadian rhythm regulation and cell cycle progression, including a novel finding of WEE1 cycling in human oral mucosa. Identified genes may serve as candidates for development of time-varying biomarkers using oral mucosal tissue collected by non-invasive brush samples. Such biomarkers, combined with other patient characteristics (such as age, sex), merit further study, particularly in cancer patients, to develop an approach for estimating individual treatment time windows for radiotherapy, a critical challenge in the application of cancer chronotherapy. Finally, although this study was motivated by reducing radiation-induced OM, results have broader implication for cancer chronotherapy in general.
The additional files for this article can be found as follows:Supplementary Figure 1
Boxplot for RIN score by sample collection time order. Samples C01-C06 were collected at approximately 10:00, 14:00, 18:00, 22:00, 2:00 and 6:00, respectively. DOI: https://doi.org/10.5334/jcr.213.s1Supplementary Figure 2
(a) Gene-feature composition plot describing genomic content distribution of oral mucosal samples by each participant and by sample collection time order. Samples C01-C06 were collected at approximately 10:00, 14:00, 18:00, 22:00, 2:00 and 6:00, respectively. Higher CDS exon content is expected from good quality RNAseq samples. (b) Percentage of mapped gene region by participant and sample collection time order. Samples C01–C06 were collected at approximately 10:00, 14:00, 18:00, 22:00, 2:00 and 6:00, respectively. DOI: https://doi.org/10.5334/jcr.213.s2Supplementary Figure 3
Gene body coverage of housekeeping genes for each collected sample by individual. Samples C01-C06 were collected at approximately 10:00, 14:00, 18:00, 22:00, 2:00 and 6:00, respectively. DOI: https://doi.org/10.5334/jcr.213.s3Supplementary Figure 4
Scatter plots showing correlation of age with gene peak times, sleep times and DLMO time. DOI: https://doi.org/10.5334/jcr.213.s4Supplementary Figure 5
Boxplot for mapping read percentage by sample collection time order. Samples C01-C06 were collected at approximately 10:00, 14:00, 18:00, 22:00, 2:00 and 6:00, respectively. DOI: https://doi.org/10.5334/jcr.213.s5Supplementary Table 1
Results of all candidate genes in the three pathways and the 108 genes genome-wide with p-value < 0.05. DOI: https://doi.org/10.5334/jcr.213.s6Supplementary Table 2
Peak time of gene mRNA expression by sex. DOI: https://doi.org/10.5334/jcr.213.s7
We thank Marina Antoch for providing advice in the protocol development. We thank the Clinical and Translational Research Center for providing facilities and staff (especially Barbara Bauer and Arpan Dharia) for the overnight sample collection. We thank the Roswell Park Data Bank and Biorepository (Warren Davis and Annmarie Nowak) for support of the development of brush cheek cell collection protocol. We thank Danielle Abramo-Balling from the Community Recruitment Liaison, Clinical and Translational Scientific Institute, and Shawna Matthews for recruiting participants.
This study was funded by the Clinical and Translational Science Institute Pilot Studies Program RF-SUNY at Buffalo/NCATS, UL1TR001412-04 (PIs: Gu and Antoch). This work was supported by National Cancer Institute (NCI) grant P30CA016056 involving the use of Roswell Park Comprehensive Cancer Center’s Data Bank and Biorepository and Genomic Shared Resources.
The authors have no competing interests to declare.
Each author made substantial contributions to the conception or design of the work; or the acquisition, analysis, or interpretation of data for the work. Each of them contributed to the drafting the work, or revising it critically for important intellectual content.
Dallmann R, Okyar A, Lévi F. Dosing-Time Makes the Poison: Circadian Regulation and Pharmacotherapy. Trends in molecular medicine. 2016; 22(5): 430–45. DOI: https://doi.org/10.1016/j.molmed.2016.03.004
Gu F, Farrugia MK, Duncan WD, Feng Y, Hutson AD, Schlecht NF, et al. Daily Time of Radiation Treatment is Associated with Subsequent Oral Mucositis Severity during radiotherapy in Head and Neck Cancer Patients. Cancer epidemiology, biomarkers & prevention: A publication of the American Association for Cancer Research, cosponsored by the American Society of Preventive Oncology 2020. DOI: https://doi.org/10.1158/1055-9965.EPI-19-0961
Blakaj A, Bonomi M, Gamez ME, Blakaj DM. Oral mucositis in head and neck cancer: Evidence-based management and review of clinical trial data. Oral oncology. 2019; 95: 29–34. DOI: https://doi.org/10.1016/j.oraloncology.2019.05.013
Lalla RV, Bowen J, Barasch A, Elting L, Epstein J, Keefe DM, et al. MASCC/ISOO clinical practice guidelines for the management of mucositis secondary to cancer therapy. Cancer. 2014; 120(10): 1453–61. DOI: https://doi.org/10.1002/cncr.28592
Sonis ST. The pathobiology of mucositis. Nature reviews Cancer. 2004; 4(4): 277–84. DOI: https://doi.org/10.1038/nrc1318
Bjarnason GA, Mackenzie RG, Nabid A, Hodson ID, El-Sayed S, Grimard L, et al. Comparison of toxicity associated with early morning versus late afternoon radiotherapy in patients with head-and-neck cancer: A prospective randomized trial of the National Cancer Institute of Canada Clinical Trials Group (HN3). International journal of radiation oncology, biology, physics. 2009; 73(1): 166–72. DOI: https://doi.org/10.1016/j.ijrobp.2008.07.009
Goyal M, Shukla P, Gupta D, Bisht SS, Dhawan A, Gupta S, et al. Oral mucositis in morning vs. evening irradiated patients: a randomised prospective study. International journal of radiation biology. 2009; 85(6): 504–9. DOI: https://doi.org/10.1080/09553000902883802
Grechez-Cassiau A, Rayet B, Guillaumond F, Teboul M, Delaunay F. The circadian clock component BMAL1 is a critical regulator of p21WAF1/CIP1 expression and hepatocyte proliferation. The Journal of biological chemistry. 2008; 283(8): 4535–42. DOI: https://doi.org/10.1074/jbc.M705576200
Kowalska E, Ripperger JA, Hoegger DC, Bruegger P, Buch T, Birchler T, et al. NONO couples the circadian clock to the cell cycle. Proceedings of the National Academy of Sciences of the United States of America. 2013; 110(5): 1592–9. DOI: https://doi.org/10.1073/pnas.1213317110
Laranjeiro R, Tamai TK, Peyric E, Krusche P, Ott S, Whitmore D. Cyclin-dependent kinase inhibitor p20 controls circadian cell-cycle timing. Proceedings of the National Academy of Sciences of the United States of America. 2013; 110(17): 6835–40. DOI: https://doi.org/10.1073/pnas.1217912110
Matsuo T, Yamaguchi S, Mitsui S, Emi A, Shimoda F, Okamura H. Control mechanism of the circadian clock for timing of cell division in vivo. Science (New York, NY). 2003; 302(5643): 255–9. DOI: https://doi.org/10.1126/science.1086271
Sancar A, Lindsey-Boltz LA, Kang TH, Reardon JT, Lee JH, Ozturk N. Circadian clock control of the cellular response to DNA damage. FEBS letters. 2010; 584(12): 2618–25. DOI: https://doi.org/10.1016/j.febslet.2010.03.017
Sancar A, Lindsey-Boltz LA, Gaddameedhi S, Selby CP, Ye R, Chiou YY, et al. Circadian clock, cancer, and chemotherapy. Biochemistry. 2015; 54(2): 110–23. DOI: https://doi.org/10.1021/bi5007354
Benloucif S, Burgess HJ, Klerman EB, Lewy AJ, Middleton B, Murphy PJ, et al. Measuring melatonin in humans. Journal of clinical sleep medicine: JCSM: official publication of the American Academy of Sleep Medicine. 2008; 4(1): 66–9. DOI: https://doi.org/10.5664/jcsm.27083
Hughey JJ, Butte AJ. Differential Phasing between Circadian Clocks in the Brain and Peripheral Organs in Humans. Journal of biological rhythms. 2016; 31(6): 588–97. DOI: https://doi.org/10.1177/0748730416668049
Plikus MV, Vollmers C, de la Cruz D, Chaix A, Ramos R, Panda S, et al. Local circadian clock gates cell cycle progression of transient amplifying cells during regenerative hair cycling. Proceedings of the National Academy of Sciences of the United States of America. 2013; 110(23): E2106–15. DOI: https://doi.org/10.1073/pnas.1215935110
Scheving LE, Burns ER, Pauly JE. Circadian rhythms in mitotic activity and 3 H-thymidine uptake in the duodenum: Effect of isoproterenol on the mitotic rhythm. Am J Anat. 1972; 135(3): 311–7. DOI: https://doi.org/10.1002/aja.1001350302
Scheving LE, Burns ER, Pauly JE, Tsai TH. Circadian variation in cell division of the mouse alimentary tract, bone marrow and corneal epithelium. Anat Rec. 1978; 191(4): 479–86. DOI: https://doi.org/10.1002/ar.1091910407
Buchi KN, Moore JG, Hrushesky WJ, Sothern RB, Rubin NH. Circadian rhythm of cellular proliferation in the human rectal mucosa. Gastroenterology. 1991; 101(2): 410–5. DOI: https://doi.org/10.1016/0016-5085(91)90019-H
Warnakulasuriya KA, MacDonald DG. Diurnal variation in labelling index in human buccal epithelium. Arch Oral Biol. 1993; 38(12): 1107–11. DOI: https://doi.org/10.1016/0003-9969(93)90173-J
Bjarnason GA, Jordan RC, Wood PA, Li Q, Lincoln DW, Sothern RB, et al. Circadian expression of clock genes in human oral mucosa and skin: Association with specific cell-cycle phases. The American journal of pathology. 2001; 158(5): 1793–801. DOI: https://doi.org/10.1016/S0002-9440(10)64135-1
Casparis S, Borm JM, Tomic MA, Burkhardt A, Locher MC. Transepithelial Brush Biopsy – Oral CDx(R) – A Noninvasive Method for the Early Detection of Precancerous and Cancerous Lesions. Journal of clinical and diagnostic research: JCDR. 2014; 8(2): 222–6. DOI: https://doi.org/10.7860/jcdr/2014/7659.4065
Duffy JF, Dijk DJ. Getting through to circadian oscillators: Why use constant routines? Journal of biological rhythms. 2002; 17(1): 4–13. DOI: https://doi.org/10.1177/074873002129002294
Zhao J, Davé SB, Wang J, Subramanian PS. Clinical color vision testing and correlation with visual function. Am J Ophthalmol. 2015; 160(3): 547–52.e1. DOI: https://doi.org/10.1016/j.ajo.2015.06.015
Burgess HJ, Wyatt JK, Park M, Fogg LF. Home Circadian Phase Assessments with Measures of Compliance Yield Accurate Dim Light Melatonin Onsets. Sleep. 2015; 38(6): 889–97. DOI: https://doi.org/10.5665/sleep.4734
Ancoli-Israel S, Cole R, Alessi C, Chambers M, Moorcroft W, Pollak CP. The role of actigraphy in the study of sleep and circadian rhythms. Sleep. 2003; 26(3): 342–92. DOI: https://doi.org/10.1093/sleep/26.3.342
Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009; 10(3): R25. DOI: https://doi.org/10.1186/gb-2009-10-3-r25
Trapnell C, Pachter L, Salzberg SL. TopHat: Discovering splice junctions with RNA-Seq. Bioinformatics (Oxford, England). 2009; 25(9): 1105–11. DOI: https://doi.org/10.1093/bioinformatics/btp120
Wang L, Wang S, Li W. RSeQC: quality control of RNA-seq experiments. Bioinformatics (Oxford, England). 2012; 28(16): 2184–5. DOI: https://doi.org/10.1093/bioinformatics/bts356
Wu G, Anafi RC, Hughes ME, Kornacker K, Hogenesch JB. MetaCycle: an integrated R package to evaluate periodicity in large scale data. Bioinformatics (Oxford, England). 2016; 32(21): 3351–3. DOI: https://doi.org/10.1093/bioinformatics/btw405
Glynn EF, Chen J, Mushegian AR. Detecting periodic patterns in unevenly spaced gene expression time series using Lomb-Scargle periodograms. Bioinformatics (Oxford, England). 2006; 22(3): 310–6. DOI: https://doi.org/10.1093/bioinformatics/bti789
Wu G, Ruben MD, Schmidt RE, Francey LJ, Smith DF, Anafi RC, et al. Population-level rhythms in human skin with implications for circadian medicine. Proceedings of the National Academy of Sciences of the United States of America. 2018; 115(48): 12313–8. DOI: https://doi.org/10.1073/pnas.1809442115
Liberzon A, Subramanian A, Pinchback R, Thorvaldsdóttir H, Tamayo P, Mesirov JP. Molecular signatures database (MSigDB) 3.0. Bioinformatics (Oxford, England). 2011; 27(12): 1739–40. DOI: https://doi.org/10.1093/bioinformatics/btr260
Hughes ME, Abruzzi KC, Allada R, Anafi R, Arpat AB, Asher G, et al. Guidelines for Genome-Scale Analysis of Biological Rhythms. Journal of biological rhythms. 2017; 32(5): 380–93. DOI: https://doi.org/10.1177/0748730417728663
Benjamini Y, Drai D, Elmer G, Kafkafi N, Golani I. Controlling the false discovery rate in behavior genetics research. Behavioural brain research. 2001; 125(1–2): 279–84. DOI: https://doi.org/10.1016/S0166-4328(01)00297-2
Bakdash JZ, Marusich LR. Corrigendum: Repeated Measures Correlation. Frontiers in Psychology. 2019; 10. DOI: https://doi.org/10.3389/fpsyg.2019.01201
Polidarová L, Sládek M, Soták M, Pácha J, Sumová A. Hepatic, duodenal, and colonic circadian clocks differ in their persistence under conditions of constant light and in their entrainment by restricted feeding. Chronobiology international. 2011; 28(3): 204–15. DOI: https://doi.org/10.3109/07420528.2010.548615
Polidarová L, Soták M, Sládek M, Pacha J, Sumová A. Temporal gradient in the clock gene and cell-cycle checkpoint kinase Wee1 expression along the gut. Chronobiology international. 2009; 26(4): 607–20. DOI: https://doi.org/10.1080/07420520902924889
Takahashi JS, Hong HK, Ko CH, McDearmon EL. The genetics of mammalian circadian order and disorder: Implications for physiology and disease. Nature reviews Genetics. 2008; 9(10): 764–75. DOI: https://doi.org/10.1038/nrg2430
Lee C, Etchegaray JP, Cagampang FR, Loudon AS, Reppert SM. Posttranslational mechanisms regulate the mammalian circadian clock. Cell. 2001; 107(7): 855–67. DOI: https://doi.org/10.1016/S0092-8674(01)00610-9
Brown SA, Fleury-Olela F, Nagoshi E, Hauser C, Juge C, Meier CA, et al. The period length of fibroblast circadian gene expression varies widely among human individuals. PLoS Biol. 2005; 3(10): e338. DOI: https://doi.org/10.1371/journal.pbio.0030338
Matheson CJ, Backos DS, Reigan P. Targeting WEE1 Kinase in Cancer. Trends in pharmacological sciences. 2016; 37(10): 872–81. DOI: https://doi.org/10.1016/j.tips.2016.06.006
Geenen JJJ, Schellens JHM. Molecular Pathways: Targeting the Protein Kinase Wee1 in Cancer. Clinical cancer research : An official journal of the American Association for Cancer Research. 2017; 23(16): 4540–4. DOI: https://doi.org/10.1158/1078-0432.CCR-17-0520
Gery S, Komatsu N, Baldjyan L, Yu A, Koo D, Koeffler HP. The circadian gene per1 plays an important role in cell growth and DNA damage control in human cancer cells. Molecular cell. 2006; 22(3): 375–82. DOI: https://doi.org/10.1016/j.molcel.2006.03.038
Roenneberg T, Kuehnle T, Pramstaller PP, Ricken J, Havel M, Guth A, et al. A marker for the end of adolescence. Curr Biol. 2004; 14(24): R1038–9. DOI: https://doi.org/10.1016/j.cub.2004.11.039
Chen CY, Logan RW, Ma T, Lewis DA, Tseng GC, Sibille E, et al. Effects of aging on circadian patterns of gene expression in the human prefrontal cortex. Proceedings of the National Academy of Sciences of the United States of America. 2016; 113(1): 206–11. DOI: https://doi.org/10.1073/pnas.1508249112