Comprehensive integrative analysis of circadian rhythms in human saliva

Data were collected from 2018 to 2023 from 21 healthy participants (52.4% female) with an average age of 31 years (25–43 years). The individual experiments consisted of 4 to 19 participants, depending on the experiment, as described below. Saliva was collected at 3–4 time points/day over 2 consecutive days for RNA extraction and core clock gene expression analysis using TimeTeller^® kits. Five participants provided saliva samples at multiple time points during the study period. Hormonal analysis (cortisol and melatonin) was conducted in six participants over 1–2 consecutive days alongside RNA sampling. PAP (Papanicolaou)-based staining of epithelial cells and leukocytes was performed in three participants concurrently with RNA sampling. Fourteen participants completed the MEQ-SA questionnaire¹² for chronotype determination (Fig. 1). For saliva volume and RNAprotect ratio optimization, four participants provided samples at a defined time point (13 h). Human, mouse, and baboon datasets (mammalian reference datasets, for details, see “Materials”) were obtained from GEO or ArrayExpress^33–35 for tissue synchrony analysis of ARNTL1 and PER2 and in-silico time series analysis to determine the minimum data points required for circadian rhythm assessment.

This study integrates diverse data types for human circadian rhythm assessment (Fig. 1), addressing technical and clinical goals. Firstly, we explored saliva’s feasibility and biological basis for the analysis of the circadian rhythm profile. We aimed to demonstrate (a) saliva’s potential as an optimal starting material, given the known tissue synchronization across clocks within the organism; (b) protocol optimization for high-quality RNA extraction from saliva; and (c) the intrinsic nature of the circadian rhythm, unaffected by cell composition variability. Importantly, saliva’s convenience for home or outpatient sampling is clinically advantageous.

Secondly, we aimed to validate our protocol’s performance and robustness for assessing circadian rhythms via clock gene RNA expression in saliva. Our goals include: (a) demonstrating assessable circadian profiles from clock gene RNA expression; (b) confirming rhythm stability over consecutive days; (c) identifying conserved correlations between circadian parameters in mouse reference profiles and our cohort; and (d) proposing strategies for sample reduction while maintaining analysis quality, beneficial for real-world applications, offering time-saving and flexibility, particularly in hospital settings.

Thirdly, we sought to correlate the obtained gene expression analysis results with established circadian rhythm assessment methods such as melatonin and cortisol hormone levels and chronotype. Clinically, we propose gene expression profiling of core clock genes as valuable for disease prevention, treatment, and diagnosis. This approach can directly link rhythmic information to metabolic networks or drug targets, facilitating personalized medicine through companion (essential for effective personalized medicine) or complementary (provides valuable information, but it is not required for the administration of the drug) diagnostics.

We focused here on specific core clock genes as they are known to steer downstream clock-controlled genes and provide a good estimate of the overall clock status, therefore providing a robust and cost-effective characterization of the expression profile of the peripheral clock³⁶. The chosen clock genes ARNTL1, NR1D1, and PER2 were shown in previous studies to exhibit a robust and well-detectable circadian rhythm in saliva and/or oral mucosa^26,28,37. Depending on the biological context, other clock genes might be meaningful. For instance, PER3, as compared to other PER isoforms, was shown to have stronger oscillation phenotypes in blood, oral mucosa, and hair follicles regarding the amplitude, inter-individual variability, and significance of rhythmicity^30,38,39. Our protocol is not limited to the clock genes shown, but can be used for other clock- and clock-controlled genes²⁶.

We previously showed the feasibility of assessing circadian rhythms in peripheral clocks in humans using saliva samples²⁶. Here, we analyzed circadian profiles of a variety of tissues in publicly available datasets of human, mice and baboon^33–35. The data showed phase synchronization of clock genes ARNTL1 and PER2 across peripheral tissues (Supplementary Fig. 1). This validates the use of peripheral tissues for circadian rhythm analysis, supporting the usage of saliva, with its practicality, ease of handling, and non-invasive collection protocols, as a suitable biological material for assessing circadian rhythms in human studies. We established a protocol for obtaining RNA in sufficient yields and quality to allow for quantification of gene expression from saliva in a previous study²⁶. To further optimize the previous protocol, we chose a study subgroup (n = 4). We used RNAprotect, a standard preservative to protect RNA from degradation, and analyzed different ratios of saliva:RNAprotect, as well as different volumes of saliva. Following RNA extraction, RNA concentration and quality/purity (as defined by A260/230 and A260/280 values) were determined. Overall, maximal yields could be obtained with a 1:1 ratio and 1.5 mL of saliva. Additionally, no significant variance in RNA quality/purity was observed across the conditions tested (Supplementary Fig. 2). Thus, we used a 1:1 ratio and 1.5 mL saliva for subsequent experiments.

Next, we investigated whether alterations in cellular composition, specifically the leukocyte/epithelial cell (L/E) ratio, might impact the circadian rhythmicity of the core-clock genes measured. For that, we collected the saliva of three participants at 3-time points (9 h, 13 h, 17 h) over 2 days, and, in parallel, assessed the cellular composition, and gene expression of ARNTL1 and PER2. The cells were stained using the PAP method, which showed good preservation of cellular properties and allowed us to distinguish between the two main cell types expected in saliva, epithelial cells, and leukocytes (Fig. 2a, b), as previously described^32,40. We could also observe bacteria from the oral flora (Fig. 2a, b) on the samples, as expected, which are unlikely to influence the RT-PCR measurements. Firstly, because no bacterial homologs to human clock genes have been found to date^41,42, secondly, because we used primers validated by the manufacturer to be specific for the human genes analyzed.

The cell composition did not show a significant time-dependency (Fig. 2c, d). Furthermore, the L/E ratio did not significantly correlate with either ARNTL1 or PER2 which indicated that the observed diurnal fluctuations of core clock genes in saliva cellular composition were indeed due to internal circadian rhythmicity of the cells rather than fluctuations in cell composition (Fig. 2e, f).

For circadian profiles, the difference between maximal and minimal expression value over time (max-min) represents an estimate of the oscillation amplitude. For the peripheral tissues of mice, we observed a correlation between max-min and mean expression level (mean over time) for all core-clock genes (Fig. 3a–c, Supplementary Fig. 3). The correlation was particularly strong for Per2, Nr1d1, and Arntl (Supplementary Fig. 3), which motivated us to measure these genes. We found a similar correlation between amplitude and mean for the gene expression in human saliva (Fig. 3d–f, full line from human data, dashed line from mouse data). The observed small reduction in slope in the human data as compared to the mouse data could be explained by the limited time period of measurements because restricting the measurements of the mouse data to a similar time period reduces the regression slope (Fig. 3a–c, green line). A correlation between max-min and mean level should imply a correlation of max-min with the maximal expression value (max), which we indeed observe in the human and mouse data (Fig. 3g–l). In this case, the match of the correlation slopes is reduced, potentially because the measurements are based on fewer data points, i.e., evaluation of the highest expression with only 4 data points per day is less robust as compared with the mean level. Similar regression lines for humans and mice suggest that we can indeed access fundamental parameters of circadian profiles, in particular, good estimates of mean expression and expression amplitude, of PER2 and ARNTL1 in our saliva-based data.

Core-clock genes interact in a transcription-translation network, and thus we investigated the correlation between the temporal mean level of those different genes. In the mouse data, we found a significant correlation over the different tissues between the expression of all pairs of genes, Per2 and Arntl, Per2 and Nr1d1, and Nr1d1 and Arntl (Fig. 3m–o). In the human data, the mean expression of PER2 and NR1D1 was correlated over different subjects and days, while PER2 and ARNTL1, as well as NR1D1 and ARNTL1, were not linearly correlated in our data set (Fig. 3p–r). The correlation in the mouse data might be individual-specific, in which case we would not expect a correlation over different human subjects, or the missing correlation for ARNTL1 might be explained by a putative stronger influence of lifestyle factors on ARNTL1 expression as compared to PER2 and NR1D1 expression. While the above correlation analysis focuses on the relative size of circadian characteristics between humans and mice, we also find that the circadian profile, when analyzing two consecutive days and taking the mean over many subjects, shows, also in absolute values, a relation to a mouse reference tissue (Supplementary Fig. 4).

Core-clock genes (ARNTL1, PER2, NR1D1) exhibit circadian variations throughout the day. We assessed the reproducibility of gene expression measurements in saliva samples across two consecutive days, using TimeTeller^® kits, for ARNTL1, PER2, and NR1D1. No significant differences were observed between Day 1 and Day 2 (Fig. 4a; Pairwise Wilcoxon test with Bonferroni correction). Some participants provided saliva samples over four consecutive days, confirming the robustness of our method, but also revealing interindividual variability in circadian parameters and profiles of ARNTL1, PER2, and NR1D1 (Fig. 4d). Despite individual differences, certain similarities emerged; for instance, ARNTL1 expression peaked after 16 h for most participants (Fig. 4b), while PER2 and NR1D1 peaked before 12 h. Additionally, saliva sampling at different times of the year showed changes in peak expression for ARNTL1, PER2, and NR1D1 in some participants (Fig. 4c). Seasonal variations, particularly in daylight duration, can directly impact individual circadian rhythms, as also previously reported⁴³. Male and female circadian rhythms differ in key aspects, such as cycle length, hormonal regulation, and sensitivity to environmental cues. Although the primary focus of this study was not to investigate sex differences—an analysis that would require age-matched participants—we did observe trends suggesting differences among male as compared to female participants in terms of mean gene expression, phase, and amplitude of circadian rhythms. While these differences were not statistically significant, likely due to the small sample size and confounding factors like age (Supplementary Fig. 5), they confirm sex-related variability in circadian regulation as previously reported³³.

For our analysis, we primarily utilized 8 saliva samples collected during the active phase of an individual over two consecutive days, from 9 h to 21 h. However, reducing the number of sampling time points for circadian rhythm assessment would enhance feasibility, streamline sampling, and reduce analysis costs. Such reduction is anticipated to improve compliance and better align sampling with the daily schedule of participants or the clinical routines, and thus it is relevant for routine measurements. To explore this, we employed a newly developed algorithm on published datasets from SW480 cells⁴⁴ and the Gtex dataset from healthy participants (transverse and sigmoid colon tissues), extrapolating circadian profiles from single-time-point samples³³, and compared circadian gene expression for ARNTL1 and PER2 gene expression. Results indicated aligned circadian expression (Supplementary Fig. 6a). Additionally, using RNA-Seq data from SW480 cells collected every 3 h over 30 h, we tested various in-silico conditions for sampling interval and frequency. Our findings suggest successful circadian rhythm assessment with at least 4 data points within a day (Supplementary Fig. 6b).

To evaluate the reproducibility of ARNTL1, PER2, and NR1D1 gene expression data, we calculated several error metrics: mean absolute error (MAE), mean squared error (MSE), root mean squared error (RMSE), and mean bias error (MBE). The MAE was 0.112, indicating an average deviation of 0.112 units between Day 1 and Day 2 expression values. The MSE was 0.021, reflecting the average squared difference, with lower values signifying better consistency. The RMSE was 0.1475, representing the typical error magnitude in the same units as the expression data. The MBE was 0.049, suggesting a slight positive bias, with Day 2 expression values slightly higher than those on Day 1. Overall, these metrics demonstrate strong reproducibility with low errors, confirming consistency in the gene expression data across days. For a more comprehensive assessment of the participants’ circadian profile, we aimed to relate the gene expression profiles determined to other methods reflecting different layers of circadian regulation, i.e., the hormones cortisol and melatonin², and chronotype assessment using the MEQ-questionnaire¹², which we measured in parallel. Spearman correlation analysis was carried out between ARNTL1, PER2, and cortisol circadian parameters (phase, amplitude), as well as between their mean levels (Fig. 5a). Cortisol mean level showed a significant negative correlation with melatonin level (r = −0.78, p = 0.04). Cortisol phase also showed a significant negative correlation with ARNTL1 phase (r = −0.75, p = 0.05) and PER2 mean expression (r = −0.75, p = 0.05; Fig. 5a). In addition, MI and MIC analyses helped in deciphering the linear and non-linear relationship between gene expression and hormonal data (Supplementary Fig. 7). MI measures the dependency between two variables, i.e., MI of 0 indicates that the variables are independent of each other and higher values indicate a stronger dependency. MIC identifies the strength of the association between the two variables, and its value closer to 1 indicates a strong association. Both MI and MIC showed that ARNTL1 expression (MI = 0.45, MIC = 0.52) and ARNTL1 phase (MI = 0.73, MIC = 0.98) had a strong relation with Cortisol phase (Supplementary Fig. 7).

We did not determine circadian profiles for melatonin as it was only measured at one timepoint per day, but found a correlation of melatonin measured at bed-time with the mean gene expression of PER2, but not with the mean gene expression of ARNTL1 (Fig. 5b). Interestingly, we found a correlation of the mean cortisol level with ARNTL1, but not with PER2 (Fig. 5b). We also found correlations between the acrophase of ARNTL1 and the acrophase of cortisol (Fig. 5c). Furthermore, both acrophases correlate with the bedtime of the subjects at the sampling day (Fig. 5c). The correlation of ARNTL1 with cortisol might hint at a stronger dependence of ARNTL1 on the particularities of daily life, such as the amount of stressor exposure as measured by the “stress hormone” cortisol. PER2 and melatonin, on the other hand, might both be correlated positively as both depend strongly on light exposure. Moreover, the circadian profile of ARNTL1, PER2, and cortisol depicted variation among different participants (Fig. 5d). Cortisol showed good rhythmicity in our participants (see Fig. 5d for representative examples).

Finally, we compared the circadian phases of ARNTL1 and PER2 with the MEQ scores of the participants. As mentioned in the introduction section, the MEQ score is seen as a rough estimate of the circadian phase¹¹. For our participants, we did not find a significant correlation between gene expression and MEQ score (Fig. 5e). The distribution of MEQ scores on our dataset (n = 14) seems to represent the same distribution as reported from a larger population study (n = 594) (Fig. 5f). This data was retrieved from Randler et al., the German version of the reduced Morningness-Eveningness Questionnaire (rMEQ)⁴⁵. Since the rMEQ correlates well with MEQ⁴⁶, and the study originates from a distribution in Germany, this data represents an appropriate reference. As expected, most participants showed an intermediate chronotype (Fig. 5f). In both populations, extreme chronotypes were not observed.

Altogether our results showed the feasibility and robustness of our protocol to measure the circadian rhythm via clock gene expression in saliva, which can be used in a health and disease context.

The data was collected between 2018 and 2023 (but not during the COVID high infection time in Germany) and consisted of 21 healthy participants with a mean age of 31 years (25–43 years) and 52.4% females. The individual experiments included 4 to 19 participants, as indicated in the results section. All participants provided saliva at 3-4 timepoints/day over 2 consecutive days for RNA extraction and gene expression analysis of core clock genes, however for participant 18, the RNA concentration was too low for gene expression analysis. Some participants provided saliva for RNA extraction more than once across the year. Six participants provided samples for hormonal (cortisol and melatonin) analysis over one or two consecutive days. Three participants gave samples for PAP-based staining of epithelial cells and leukocytes in saliva. The sampling for the hormonal and cellular measurements was performed in parallel to that for RNA extraction and gene expression analysis. In parallel, each participant filled out the Morningness-Eveningness-Questionnaire-Self Assessment (MEQ-SA) questionnaire for chronotype determination¹². For the optimization of saliva volume and saliva/RNAprotect ratio for RNA extraction a total of four participants provided samples at a defined timepoint (13 h) to avoid possible time-of-day dependent fluctuations of saliva composition and RNA amount in saliva.

All procedures involving human subjects were approved by the Charité—Universitätsmedizin Berlin Ethics Committee (EA2/188/18; EA4/040/20). The participants were enrolled on a voluntary basis and signed the informed consent approving the analysis of their data in the scope of the study.

Unstimulated saliva samples were collected using TimeTeller^® Kits at 3-4 timepoints (9 h, 13 h, 17 h, 21 h) for 2 consecutive days. One hour before sample collection participants refrained from eating, drinking, smoking, and oral hygiene. Saliva was collected in 3 mL collection tubes with RNAprotect Cell reagent (Cat No. 76526, Qiagen), the tube was shaken vigorously to mix and stored at 4 °C for up to 1 month or at -20 °C for long-term storage. For ratio optimization, the ratios saliva:RNAprotect 1:1 @ 1.5 mL saliva, 1:2 @ 1 mL saliva, 2:1 @ 1 mL saliva, and 1:2 @ 0.5 mL saliva were tested. For all other experiments, the ratio was 1:1 (1.5 mL saliva and 1.5 mL RNAprotect).

The participants who gave consent for the collection of their hormonal levels were provided with a commercially available kit for Cortisol and Melatonin analysis from Cerascreen®. For cortisol measurements, saliva was sampled directly after and 0.5 h, 1 h, 2 h, 5 h, 8 h and 12 h upon awakening. For melatonin quantification, saliva was sampled 30 min before going to bed.

RNA was extracted as previously described using TimeTeller® Kits according to the manufacturer indications^9,26. Gene expression fold-change was calculated using the ΔΔCT method normalizing CT values first to GAPDH and then to mean expression value of the respective target gene. The procedures have been described elsewhere in detail²⁶.

50 µL saliva was pipetted on a glass microscopy slide and immediately fixed using RotiFix spray (Cat No. CL85.1, Roth, Karlsruhe, Germany). After drying the fixing spray, the slides were further processed.

For PAP staining the PAP rapid stain kit (Cat No. 9195.1, Roth, Karlsruhe, Germany) was used according to the instructions of the manufacturer. Briefly, unstimulated saliva was collected for 2 consecutive days at 9 h, 13 h, and 17 h 1.5 mL saliva was mixed with 1.5 mL RNA protection for gene expression analysis as described above. 50 µL saliva was spread on a glass microscopy slide and fixed using RotiFix Spray. After drying, the slides were incubated in a descending ethanol series (96%, 80%, 70%; 1 h each). Afterward, slides were incubated in solution 1 (modified hematoxylin solution; 3 min), blued in tap water for 3 min, and incubated in an ascending ethanol series (70%, 80%, 90%; 30 s each). Slides were then incubated in solution 2 (modified polychrome solution), rinsed in 96% ethanol, and dried in 96% ethanol (2 × 30 s) and 100% ethanol (3 × 30 s). Roticlear (Cat. No A538.1, Roth, Karlsruhe, Germany) was used as Intermedium (2 × 1 min, 1 × 3 min), and slides were mounted in RotiMount (Cat. No HP68.1, Roth, Karlsruhe, Germany) and coversliped. After hardening of the mounting medium slides were stored at 4 °C.

Gene expression data was processed using the multivariate harmonic regression model E_gc = μg + a g cos (φ c) + bg sin (φc) + ε, where “c” represents the condition, “g” represents the gene, and “ε” denotes the Gaussian noise represent random variability in the gene expression measurements as also previously used³³. The model fits a periodic function (using cosine and sine terms) to the gene expression data, accounting for a 24-h cycle (a fixed period of 24 h used in our analysis). Likelihood Ratio Test is then used to compare the two models (with and without periodic terms) to estimate amplitude and phase. The same procedure was used for the RNA-seq and microarray datasets analyzed. Phase is a circular variable, with phases before and after midnight close in time.

For cortisol, we used a harmonic regression with soft L1 cost, fitting μ + a cos (φ) + b sin (φ) to the cortisol measures, using the Python function scipy.optimize.least_squares(error, x0 = [a = −0.31084448, b = 2.33838041, μ = 6.33268443], loss = ‘soft_l1’, f_scale = 0.1) with error the absolute difference between experimental data and fit.

For the mouse microarray dataset (GSE54650↗) from Zhang et al.³⁴, publicly available data was retrieved from the Gene Expression Omnibus (GEO) database and pre-processed using the Robust Multichip Average (RMA) pipeline as described⁸¹. The baboon dataset (GSE98965↗) was downloaded from the European Nucleotide Archive (ENA), and SW480 and SW620 RNAseq data is available under the ArrayExpress database under the accession number E-MTAB-7779⁴⁴. RNA-sequencing reads from both baboon and cell line data were aligned to reference genomes, using STAR aligner with default parameters⁸¹. Gene summarized counts were then normalized using the trimmed mean of M-values (TMM) method, and counts were log2-transformed to the counts per million (CPM) scale. A low gene expression filter of at least 0.5 CPM on average over the entire time series (per tissue) was applied before downstream analysis.

For a direct comparison with our human data, we loaded the series matrix file of the mouse microarray dataset (GSE54650↗) from Zhang et al.³³. The microarray expression values of Arntl, Per2, and Nr1d1 were normalized at each timepoint by the expression value of Gapdh, the housekeeping gene used for the qPCR analysis, which allows us to compare the absolute values of microarray data and qPCR data.

For the correlation analysis of the mice data (Fig. 3, blue), we evaluated the Gapdh-normalized data by calculating the amplitude as the difference between minimal and maximal expression value, the maximum expression as the maximal value, and the mean level as the mean overall values over the two days for each peripheral tissue (aorta, adrenal gland, brown fat (anterior dorsum adipose), heart, kidney, liver, lung, skeletal muscle (gastrocnemius), and white fat (epididymal adipose)). For the analysis restricted to the human sampling interval (Fig. 3, green), we evaluated the same for the data measured at 4h-14h and 26h-38h since wakeup (see Supplementary Fig. 4). For the correlation analysis of the human data (Fig. 3, violet), we evaluated the gene expression profiles for each subject and day separately, assuming statistically independent measurements for each day for the linear regression analysis. For each set of up to 4 sampling values per day, we calculated the amplitude as the difference between minimal and maximal expression values, the maximum expression as the maximal value, and the mean level as the mean overall values associated with this day and subject.

For Supplementary Fig. 4, we plotted data from different subjects for two consecutive days. While some subjects contributed multiple pairs of days, the calculation of the standard error of the mean (SEM) assumes statistical independence between each pair of days. We excluded data with two or less datapoints over the two days. For Fig. 5, hormonal values for cortisol that were reported as “<2” were plotted for representative purposes as 1, and one melatonin value reported as “>100” was omitted from the analysis.

Pairwise Spearman correlation was carried out between ARNTL1, PER2, and Cortisol circadian properties (amplitude, phase), melatonin and cortisol levels, and mean expression (2^dcq) using R. For Fig. 3 and Fig. 5b, c, we calculated the regression line with scipy.stats.linregress() with standard parameters. Reported are the two-sided p-values of this function. To analyze the linear and non-linear relationship between gene expression data (ARNTL1, PER2) and hormonal data, we carried out mutual information (MI) and maximal information coefficient (MIC) analysis using “infotheo” and “minerva” packages in R.

Wilcoxon paired test followed by Bonferroni’s correction was carried out to evaluate the data reproducibility between the two days of saliva sampling. To assess the robustness of our method, statistical error metrics such as mean absolute error (MAE), mean squared error (MSE), root mean squared error (RMSE), and mean bias error (MBE) were carried out. MAE measures the average magnitude of the errors between Day 1 and Day 2, without considering their direction (i.e., whether the errors are positive or negative). MSE is the average of the squared differences between Day 1 and Day 2. It gives more weight to larger differences because the errors are squared before averaging. RMSE is the square root of the MSE and gives an error value in the same units as the gene expression data. MBE shows the direction of the error, indicating whether, on average, the values from Day 2 are higher or lower than those from Day 1. The gene expression data from Day 1 was used as a “predicted” value and from Day 2 as an “observed” value to calculate the error metrics.

Data is provided within the manuscript, or the supplementary information, or can be obtained upon request from the corresponding author, and following the rules of the ethics procedures of our university.