Diurnal oscillations in human salivary microRNA and microbial transcription: Implications for human health and disease

The microbiome plays a vital role in human health and disease. Interaction between human hosts and the microbiome occurs through a number of mechanisms, including transcriptomic regulation by microRNA (miRNA). In animal models, circadian variations in miRNA and microbiome elements have been described, but patterns of co-expression and potential diurnal interaction in humans have not. We investigated daily oscillations in salivary miRNA and microbial RNA to explore relationships between these components of the gut-brain-axis and their implications in human health. Nine subjects provided 120 saliva samples at designated times, on repeated days. Samples were divided into three sets for exploration and cross-validation. Identification and quantification of host miRNA and microbial RNA was performed using next generation sequencing. Three stages of statistical analyses were used to identify circadian oscillators: 1) a two-way analysis of variance in the first two sample sets identified host miRNAs and microbial RNAs whose abundance varied with collection time (but not day); 2) multivariate modeling identified subsets of these miRNAs and microbial RNAs strongly-associated with collection time, and evaluated their predictive ability in an independent hold-out sample set; 3) regulation of circadian miRNAs and microbial RNAs was explored in data from autistic children with disordered sleep (n = 77), relative to autistic peers with typical sleep (n = 63). Eleven miRNAs and 11 microbial RNAs demonstrated consistent diurnal oscillation across sample sets and accurately predicted collection time in the hold-out set. Associations among five circadian miRNAs and four circadian microbial RNAs were observed. We termed the 11 miRNAs CircaMiRs. These CircaMiRs had 1,127 predicted gene targets, with enrichment for both circadian gene targets and metabolic signaling processes. Four CircaMiRs had “altered” expression patterns among children with disordered sleep. Thus, novel and correlated circadian oscillations in human miRNA and microbial RNA exist and may have distinct implications in human health and disease.


Introduction
The proper regulation of sleep in humans is critical for normal mental and physical health. Most major organ systems exhibit fluctuations in their functional state related to sleep-wake cycles or circadian rhythm [1][2][3]. Disturbances in sleep or disruption of circadian rhythm are a common problem in many chronic brain disorders, including autism, depression, Parkinson's, and Alzheimer's. These symptoms have a negative impact on activities of daily living [3].
During sleep-wake cycles there are numerous molecular, cellular, and physiological changes that occur. Many of these changes are driven by circadian regulatory genes, such as CLOCK and BMAL [4]. These, in turn, cause a vast array of changes in the expression of physiologically significant genes, proteins, and hormones, influencing nearly every body system. However, apart from light-dark cycles, the factors that influence expression of circadian rhythm are not fully understood.
MiRNAs are found in nearly all body cells, tissues, and biofluids [10,14]. Because miRNAs regulate the majority of human genes, a considerable number of circadian genes are now thought to be directly under their influence, including CLOCK and BMAL, among others [15]. MiRNAs that circulate throughout the body in extracellular fluids are also resistant to enzymatic degradation [16], and thus may act as critical components of a molecular endocrine system [17]. Indeed, there are now considerable data implicating miRNAs in the control of various endocrine and metabolic tissues, such as the pineal and pituitary glands [18], the hypothalamus, and the gastrointestinal (GI) tract. Furthermore, disruption of circadian regulation by miRNAs can lead to significant pathology [19].
Notably, the activities of miRNAs in the gut appear to extend beyond the regulation of host gene expression, and include a strong relationship with the resident bacteria of the microbiome [20,21]. Within the GI system, the microbiome contributes to energy harvesting by generating numerous metabolites and intermediates that influence the function of other organ systems, including the brain and endocrine organs [22]. Recent evidence also indicates that there are circadian changes in the gut microbiome [23]. Thus, cross-talk between host miR-NAs and the GI microbiome may work in concert to influence temporal changes in gene expression that drive host behavior and disease.
To our knowledge, only one prior study has demonstrated diurnal variations for a select number of cell free microRNAs in human plasma, using quantitative RT-PCR [24]. However, no prior studies have harnessed next-generation sequencing to investigate diurnal variations for the entire micro-transcriptome, or explored these diurnal patterns in the GI tract parallel to the microbiome. We hypothesized that 1) a saliva-based collection method would identify host miRNA and microbial RNA elements with consistent and parallel circadian oscillations; 2) these RNA elements would target functionally-relevant biologic pathways related to host immunity, circadian rhythm, and metabolism; and 3) a subset of circadian miRNAs would demonstrate "altered" expression in a cohort of children with disordered sleep patterns.

Salivary miRNA analysis
An overview of the sample sets and analyses is provided (Fig 1). Sample set 1 contained 24 saliva samples collected at 2 time-points (~9 AM, 9 PM) on 3 days from 4 participants. There were a total of 98 miRNAs in set 1 with a significant effect of collection time (FDR < 0.01) and no effect of day of collection (FDR > 0.05). Sample set 2 contained 48 samples collected at 4 time-points (~9 AM, 1:30 PM, 5:30 PM, 9 PM) on 4 days from 3 participants. There were a total of 123 miRNAs in set 2 that showed a significant effect of collection time and no effect of day. Levels of 61 miRNAs were similarly affected by time of collection in both sample sets and were defined as putative CircaMiRs (S1 Table).
Hierarchical (heat map) clustering using salivary concentrations of the 61 CircaMiRs was performed for sample set 1 (Fig 2A) and sample set 2 ( Fig 2B). In both sample sets, the majority of CircaMiRs (n = 49; 80%) demonstrated lower levels in the morning and higher levels in the evening. Examination of the 61 miRNAs across four time points (sample set 2) revealed only a single oscillation (i.e. a single daily peak) between 9AM and 9PM. These daily oscillations were consistent across days of collection and across participants, as reflected by the lack of significant day effects in the two -way ANOVA (S1 Table).
From the 61 CircaMiR candidates, 11 miRNAs were identified as robust multivariate predictors of collection time through a feature selection algorithm using a linear regression analysis. The regression model accurately predicted collection time in all 3 sample sets, with Multiple R values ranging from 0.805-0.956 and Adjusted R 2 values ranging from 0.54-0.833 (Table 1). Notably, the multivariate model performed best when applied to samples collected during a wakeful state (9 AM-12 AM) and model performance significantly improved in sample set 3 when 4 AM samples were excluded (Adjusted R 2 = 0.880 vs 0.794, Table 1, upper). This improvement was due to non-linear trends in the expression data during the overnight period (a circadian oscillation of high values back to low values and vice versa). In fact, the predictive utility of the linear regression model (R 2 = 0.79; Fig 3A) was even found to be inferior to a non-linear regression model that used the sine-transformed average miRNA values for just one of the 61 CircaMiRs in the set 3 samples (R 2 = 0.93; Fig 3B). Interestingly, further inspection of the alpha (intercept) and beta (slope) coefficient terms across the independent sample set regressions indicated a very high degree of internal consistency in these models (Table 2), with highly significant correlations present between all sets of model term comparisons except sample set 1 and sample set 3 with the 4 AM samples included.

Salivary microbiome analysis
Sample set 1 contained a total of 75 microbial RNAs with a significant effect of collection time (FDR < 0.01) and no effect of day of collection (FDR>0.05). Sample set 2 contained a total of 32 microbial RNAs with a significant effect of collection time and no effect of day of collection. Eleven microbial RNAs with diurnal oscillations in sample sets 1 and 2 overlapped ( Table 3). The 11 RNAs from these 11 distinct microbial species were defined as putative CircaMicrobes, and examined for their ability to predict collection time in sample set 3.
A multivariate linear regression model utilizing the 11 microbial RNAs was also able to accurately predict collection time in all 3 sample sets, with Multiple R values ranging from 0.770-0.927 and Adjusted R 2 values ranging from 0.468-0.732 (Table 1). As with the miRNA model, a non-linear relationship between the time of collection and microbial RNA concentrations in sample set 3 reduced the overall accuracy of the microbial model across the full 24 hour time cycle compared to when the 4 am samples were removed from analysis (Adjusted R 2 = 0.468 vs 0.624, Table 1), which yielded results comparable to those seen in sample sets 1 and 2. Likewise, inspection of the alpha (intercept) and beta (slope) coefficient terms across the independent sample set regressions again indicated a very high degree of internal consistency in these models with highly significant correlations present between all sets of model term comparisons (Table 2).

Relationship between CircaMiRs and CircaMicrobes
Relationships between levels of the 11 CircaMiRs and the 11 microbes with oscillating transcriptional activity were assessed across all 120 samples from sample sets 1, 2, and 3 using a Pearson's correlation analysis. With the exception of one CircaMiR (miR-200b-3p) and one CircaMicrobe (Macrococcus caseolyticus), the CircaMiRs and CircaMicrobes were generally

CircaMiR target genes
Functional analysis of the 11 CircaMiRs in DIANA miRPath revealed 1265 high confidence (p<0.05, Micro-T threshold ! 0.95) mRNA targets with enrichment for 22 KEGG pathways (Table 4). Notably, 11/22 KEGG pathway targets were involved in cell signaling. Interestingly, circadian rhythm was not among the KEGG pathways targeted by the 11 CircaMiRs according to this analysis. However, of the 30 human mRNAs in the circadian rhythm KEGG pathway (hsa04710), four (13%; Csnk1e, Rora, BHLHE40, and Prkaa2) were targeted by the 11 Circa-MiRs. To more closely examine the potential relationship of CircaMiRs and circadian function, we expanded the analysis to the initial 61 CircaMiRs and used IPA software (which included additional circadian mRNA targets). The results revealed a significant overlap in Circadian Rhythm Signaling targets (13/34 mRNAs, 38%, p = 2.2e-38) based on moderate-to-high probability predicted interactions, or experimentally-observed interactions. A complete list of the 28 mRNA transcript isoforms encompassing the 13 mRNAs and their 37 CircaMiR interactors is provided (S2 Table). The MiRpath target mapping tool also failed to detect enrichment of KEGG pathways involved in immune function or bacterial regulation among the 11 CircaMiR targets (or the 5 CircaMiRs with microbial associations in Fig 4). However, several of the CircaMiRs that mapped to circadian genes were found to target mRNAs that were clearly involved in immune function (S2 Table). Subsequent interrogation of the protein-protein interaction network for all 1127 unique mRNA targets of the 11 most robust CircaMiRs using STRING software, revealed 3794 edges (interactions) with a clustering coefficient of 0.32. This exceeds the number of protein-protein interactions expected by chance alone (p = 1.0E-16), and implies interrelatedness of CircaMiR targets. Among the expected protein targets, 471 were involved in regulation of metabolic process (GO:0019222; FDR = 8.5E-23), 413 were involved in regulation of macromolecule metabolic process (GO:0060255; FDR = 4.9E-22), and 425 were involved in regulation of cellular metabolic process (GO:0031323; FDR = 8.9E-22; S3 Table).

Metabolic targets of the oral microbiome
RNA expression patterns of oral microbes from the 9 participants in sample sets 1, 2, and 3 were examined for evidence of diurnal variations in metabolic and functional clusters across four time periods: 7-9 AM, 10 AM-2 PM, 3-6 PM, and 7-10 PM. Among the 202 functional clusters targeted by microbial RNAs, 22 pathways demonstrated nominal (p<0.05) differences in representation across the four time periods (Fig 5A). Four of these functional pathways (nucleotide sugar biosynthesis, galactose; replication recombination and repair; sphingolipid metabolism; and purine metabolism) survived multiple testing corrections (FDR 0.15). Among the 22 functional pathways with nominal changes, a cluster of seven pathways was upregulated mid-day (10 AM-2 PM), while 10 pathways demonstrated diurnal peaks in the morning (7-9AM) and evening (7-10 PM). Visualization of functional pathway expression differences in a partial least squared discriminant analysis resulted in partial separation of the four time periods, while accounting for 20.6% of the variance in COG/KEGG data in two dimensions ( Fig 5B).
Based on the ability of the 11 CircaMiRs to predict time of collection in 11 typically developing, healthy children (and adults) in sample sets 1, 2 and 3, we also used a multivariate regression model examining their ability to predict time of collection in the 63 children with ASD and a normal sleep pattern, and the 77 children with ASD and comorbid disordered sleep. As we had seen in sample sets 1, 2 and 3, these 11 CircaMiRs yielded a significant regression (R 2 = 0.41, F 1,11 = 3.19, p <0.0023) that accurately predicted the time of collection with a mean absolute error of 15.25% (Fig 7). Inspection of the multivariate regression coefficients   Temporal oscillation in the human oral microtranscriptome and T scores indicated that no individual miRNA was significant in the presence of the others, although three showed strong trends (miR-629-5p, miR-22-5p, and miR-128-3p) ( Table 6). In contrast to the significant regression findings for the ASD children without sleep disorder (n = 63), the regression results for the ASD children with sleep disorders (n = 77) using the 11 CircaMiRs did not yield a significant result (R 2 = 0.20, F 1,11 = 1.46, p >0.167).

Relationships of CircaMiRs and CircaMicrobes with daily activities
Pearson's correlation analysis was used to explore relationships between oscillating salivary RNAs and three daily routines (sleep, tooth brushing, and eating) in sample set 3. Levels of 3 CircaMiRs and 5 CircaMicrobes were significantly (FDR<0.05) associated with time since last sleep (measured in hours; Table 7). Levels of all five CircaMicrobes were inversely associated with time since sleep, while 2/3 CircaMiRs were positively correlated with time since last sleep. There were 4 CiraMiRs and 5 CircaMicrobes associated with time since last tooth brushing. Levels of all five CircaMicrobes were inversely associated with time since last tooth brushing, while 3/4 CircaMiRs were positively associated with time since last tooth brushing. MiR-200b-3p was the single CircaMiR inversely associated with both sleep and tooth brushing. Notably, expression patterns of miR-200b-3p were also hierarchically clustered with CircaMicrobe expression (Fig 4). Two Prevotella and two Fusibacterium CircaMicrobes were associated with both sleep and tooth brushing. There was only one CircaMicrobe positively associated with time since last meal (and 0 CircaMiRs).

Discussion
In the present study, 61 total human miRNAs (CircaMiRs) and 11 total microbes (CircaMicrobes) displayed consistent diurnal oscillations in saliva samples obtained from 9 different children and adults collected across multiple days and times. From these, 11 miRNAs and 11 microbes were capable of accurately and reliably predicting time of sample collection. Diurnal levels of five CircaMiRs and four CircaMicrobes were strongly associated with one another. Functional analyses of the circadian RNA components displayed enrichment for numerous signaling mechanisms, particularly metabolic pathways. However, CircaMiR and CircaMicrobe levels were more strongly associated with sleep routines than with eating routines. This may explain, partly, why levels of four CircaMiRs were "altered" in autistic children with disordered sleep, relative to autistic peers without a sleep disorder and why a portion of these Circa-MiRs target circadian mRNAs. Six of the oscillating miRNAs identified in this study (miR-15b-3p, miR-24-3p, miR-106b-3p, miR-140-3p, miR-150-5p, miR-203a-3p) were among the 26 plasma miRNAs previously found to have diurnal variations in peripheral blood samples from healthy individuals [24]. These overlapping miRNAs from two distinct biofluids may represent either primary regulatory elements or primary readouts of circadian rhythmicity. This premise is supported by the fact that levels of both miR-24-3p and miR-203a-3p are disrupted in the cohort of autistic children with disordered sleep patterns. In addition, we found suggestive evidence that miR-203a-3p was associated with sleep initiation difficulties (R = 0.20; p = 0.034). Another CircaMiR, miR-142-5p, targets the clock gene RORA. Notably, miR-142-5p also displays correlated diurnal expression with its mRNA targets NAMPT (whose gene product modulates circadian clock function by releasing the CLOCK/ARNTL/BMAL heterodimer [25]) and GRIN2B (whose gene product encodes the NR2B subunit of the NDMA receptor essential to MAPK signaling in the suprachiasmatic nucleus and CaMK II signaling in the hippocampus [26]). Notably, a well-described developmental switch from NR2A to NR2B subunit expression is considered a hallmark of synaptic maturation that promotes memory formation, and elevation in miR-142-5p (which would suppress NR2B expression) is associated with amyloid beta pathology in postmortem brain samples of subjects with Alzheimer's disease (AD) [27]. The importance of this finding is highlighted by the fact that AD is associated with significant circadian pathology (e.g. "sundowning") and that miR-142-5p restores normal synapse formation and maturation (as measured by PSD95 expression) in differentiated neural cultures [28]. Temporal oscillation in the human oral microtranscriptome Such a mechanism might even contribute to the recently described circadian oscillation in synaptic spine number that has been described across different species, especially dendritic spines on inhibitory neurons in multiple brain regions [29][30][31]. Circadian miRNAs found in both plasma and saliva may also direct diurnal physiologic processes common to both peripheral biofluids. Indeed, mapping of KEGG pathway targets for the six overlapping miRNAs reveals enrichment for broad signaling mechanisms such as Wnt Signaling, Rap1 signaling, and Endocrine factor-regulated calcium reabsorption. This is consistent with functional analysis of the 11 CircaMiRs, which also display enrichment for Rap1 and other broad signaling processes (Ras, ErbB, PI3K-Akt, mTOR, MAPK). Salivary Cir-caMiRs also demonstrate target enrichment for endocrine factors (estrogen and prolactin signaling), which regulate peripheral physiologic processes in a circadian manner [32].
Unlike oscillating miRNAs in plasma, the CircaMiRs and CircaMicrobes in saliva appear uniquely geared toward metabolic functions. CircaMiR targets display enrichment for lysine degradation, choline metabolism, and insulin signaling ( Table 3). The protein products of these mRNA targets also exhibit enhanced biologic interaction in metabolism at the cellular and macromolecule levels (S3 Table). Specifically, interactions between miR-130b-3p/ATXN2, and miR-142-5p/NAMPT (Table 4) may play important roles in regulation of host metabolism, given that loss of function mutations in both ATXN2 and NAMPT are associated with obesity and diabetes mellitus [33,34].
Oscillating RNA expression within the oral microbiome also shows relationships with diurnal metabolism. Microbial RNAs appear to target KEGG and COG pathways in a diurnal manner, by up-regulating RNAs involved in terpenoid biosynthesis, gluconeogenesis, pentose phosphate pathways, and carbon fixation during the morning and afternoon time periods. In comparison, pathways related to cell replication, nucleotide biosynthesis, and purine metabolism demonstrate both morning and evening peaks. Thus, as a whole, the oral microbiome may have evolved energy utilization patterns that capitalize on the timing of host meals to extract biosynthetic materials and allow for night time replication. Interestingly, however, levels of the 11 CircaMicrobes do not appear to correlate with time since last meal. Thus, these 11 individual entities may serve a more commensal function whose metabolic activities aid host circadian rhythms. Indirect evidence for this may be found in the circadian rhythm of terpenoid biosynthesis (Fig 5A), a diverse class of hydrocarbons present in plant-based cannabinoids, or anti-inflammatory curcuminoids that play an essential role in steroid production [35]. Given the well-established rhythmicity of steroid production, this is one mechanism by which the microbiome may contribute to host circadian biology [36].
Further evidence for a synergistic relationship between CircaMicrobes and human hosts is found in the strong associations among CircaMiR and CircaMicrobe expression (Fig 4). It is somewhat surprising that CircaMiRs have few immune, or antimicrobial targets. However, this may be because the circadian components of the oral microbiome serve a commensal function. The majority of CircaMicrobes are not known to play pathogenic roles in human hosts. Of the 11 CircaMicrobes, only three are distinct human pathogens (Haemophilus parainfluenza T3T1, Haemophilus, and Campylobacter hominis ATTC BAA-381) and none of these three are associated with CircaMiR levels. Instead CircaMiRs may interact with the oral microbiome to coordinate metabolic patterns, or production of essential amino acids. Perhaps metabolic activity by the oral microbiome leads to changes in host miRNAs that regulate downstream physiologic pathways.
To our knowledge this is only the second study to report consistent circadian rhythmicity of peripheral miRNA expression in humans, and the first to do so in saliva. The importance of this finding is underscored by the vast number of publications seeking to use peripheral miRNA expression as a biological marker of human disease [37], a venture that could be greatly confounded by failure to control for time of collection. For example, among seven studies describing peripheral miRNA expression (saliva, serum, lymphoblasts [38]) in patients with ASD, none have reported, or controlled for time of RNA collection. These studies have reported a combined 139 ASD-related miRNAs. Notably, 10 of these (7%) overlap with the 61 CircaMiRs identified herein, which could represent confounded results. Future studies may be able to utilize CircaMiR levels to control for circadian variation in miRNA expression or accurately identify time of collection among biofluid samples.
The current study also adds to the growing body of literature that suggests miRNAs may serve as a communication mechanism between the gut microbiome and human hosts [39]. Specifically, these results show how miRNA-microbiome cross-talk may occur in a circadian manner. Given the diurnal rhythmicity of human metabolism, this finding has implications in human health and disease. For instance, daily fluctuations in host-microbiome interaction may inform risk for obesity, or insulin resistance (an enriched KEGG target of the 11 Circa-MiRs). Alternatively, disruptions in miRNA-microbiome networks may unsettle the gutbrain-axis, a concept implicated in diseases such as Parkinson's [40] and ASD [41] (both of which are associated with disordered sleep).
There are several notable limitations to this study that must be considered in interpreting these findings. Detailed daily activity logs were available only for the participants in sample set 3. The remaining participants reported no medical co-morbidities (including disordered sleep), though timing of sleep initiation and cessation were not recorded. Such information, when recorded alongside physiologic measurements such as sleep architecture, or melatonin flux could be extremely informative when interpreting RNA results in future studies. Nonetheless, the RNA expression patterns from participants in sample sets 1 and 2 were sufficient to accurately predict collection time in a third independent sample set with documented sleep wake cycles.
Notably, predictive performance in sample set 3 was somewhat impaired for the subset of samples obtained at 4 AM. This may have resulted because the sinusoidal model created from samples collected between 8 AM and 8 PM could not fully account for the overnight rhythmicity that occurs in a sleep state. There may also be microbial variability introduced by differences in participant breathing patterns (e.g. open-mouthed versus nasal breathing) or fasting during sleep. Certainly, a more controlled study which tightly dictated wake time, sleep initiation time, diet, dental hygiene, and other factors could account for time of collection with greater precision. However, the current results demonstrate that even in the face of typical variability among daily routines, these 11 miRNAs and 11 microbial RNAs are remarkably accurate predictors of time of saliva collection in four different cohorts of human subjects.
The accuracy of these results may even be underestimated given the broad age range (3-55) of participants in sample sets 1, 2, and 3. The CircaMiR and CircaMicrobe candidates were generated from 2 cohorts of children and validated in a cohort of teens and adults. This is despite the fact that teens are known to have altered circadian rhythm compared with pre-teen peers and adults [42]. Circadian RNAs from sample sets 1-3 also demonstrated significant relationships with collection time in a large cohort of children with ASD. Thus, the age and developmental diversity of these sample sets may be viewed as a confounding variable, but it likely enhances the veracity of these results.
Finally, it should be noted that the RNAseq approach used to identify oral microbes and estimate transcriptional activity of individual taxons differs from the typical 16S approach used to measure microbial abundance. Thus, these results should not be interpreted as diurnal fluctuations in the quantity of the oral microbiome, but rather as circadian variation in salivary microbial activity. RNAseq and 16S measures are complementary (though not equivocal) and could potentially add to the interpretive value of this approach in future studies. Such studies might also utilize animal models to explore the cellular origins of salivary CircaMiRs, and investigate the mechanisms regulating CircaMiR production, transport, and degradation. Manipulating the gut microbiome in this setting may also provide insights into microbial-miRNA communication.
Parallel circadian oscillation in host and microbial RNA represents an important consideration for studies analyzing epi-transcriptomic or metagenomic mechanisms in human health and disease. Circadian rhythm disturbances are a common problem in disorders of the central nervous system (e.g. Parkinson's, Alzheimer's, autism, depression, concussion [43]). Hence, studies of peripheral miRNA expression in these conditions might consider how diurnal miRNA expression patterns are shifted, rather than simply focusing on average miRNA levels at a single collection point in comparison with a control cohort. Monitoring levels of these factors in biofluids like saliva could have diagnostic potential in diseases with altered circadian rhythm and may one day provide a basis for targeted miRNA therapy of circadian disruptions.

Subject assessment
This study was approved by the Institutional Review Board for the Protection of Human Subjects (IRB) at SUNY Upstate Medical University. Informed written consent was obtained for nine healthy human volunteers, and verbal assent was provided by all participating children.

Study design
A prospective cohort design employing high throughput RNA sequencing was used to examine salivary RNAs (human and microbial) for daily oscillations in concentration (Fig 1). Nine healthy participants (3-55 years of age) were divided into three groups, and provided multiple saliva samples across a unique multi-day timeline (described below). Overlapping circadian RNA candidates from the first two independent sample sets were validated in a third sample set. Human miRNAs and microbial RNAs with confirmed diurnal variation were examined for associations in expression levels. Relationships between oscillating miRNAs and coding mRNA targets were also explored. Finally, the circadian RNA components were interrogated for functional relevance to human health and disease with the following three steps: 1) mRNA target networks for human miRNAs were identified in DIANA miRPath and Ingenuity Pathway Analyst (IPA, Qiagen), while metabolic pathways targeted by microbial RNAs were defined with MicrobiomAnalyst; 2) oscillating RNAs were retrospectively interrogated in a cohort of 140 children with autism spectrum disorder (ASD) with comorbid (n = 77), or absent (n = 63) sleep disturbance; and 3) the relationship of diurnal salivary RNAs with daily activities (tooth brushing, sleep, and eating) was assessed through Pearson correlation testing.
This study examined human miRNA and microbial RNA in saliva, because this biofluid provides on-demand access to repeated sampling of the GI tract at its sole point of entry, and represents a major site of host-environment interaction. Furthermore, studies of salivary miRNA in human patients have previously shown connections with brain-related dysfunction and potential relationships with time of collection [44,45].

Participants
Participants included nine healthy volunteers, taking no daily medications, with no history of hospitalization, surgery, or sleep disorder. None of the participants had active dental caries. The nine participants were 3-55 years of age, 55% male, and 100% Caucasian. Participants provided saliva samples at various times of day on repeated days in four different sets of samples: Sample Set 1: Morning and evening samples (n = 24) collected at approximately 9 AM and 9 PM on days 1, 3, and 7 for 4 children (two male, two female; average age 7.5 yrs); Sample Set 2: Early morning, early afternoon, late afternoon, and early evening samples (n = 48) collected at approximately 9 AM, 1:30 PM, 5:30 PM, and 9 PM on days 1, 5, 10 and 15 for three female children (average age 5.1 yrs), of whom two were part of Sample Set 1; Sample Set 3: 12 samples collected at various times (ranging from 4 AM to midnight) on days 1 and 2 on two male children (average age 16.0 yrs) and their male and female parents (average age 51.5 yrs). Notably, detailed data regarding time of sleep, meals, and tooth brushing was collected for participants in Sample set 3.
Sample Set 4: Functional analysis of circadian RNAs was performed through retrospective analysis of data from an additional cohort of 140 children with ASD and comorbid sleep disturbance (n = 77), or normal sleep (n = 63). Salivary RNA was collected from these 140 children (2-6 years of age) at a single time-point, between 8 AM and 8 PM in a non-fasting state. ASD was confirmed through physician diagnosis, using the Diagnostic and Statistics Manual of the American Psychiatric Association, 5 th Edition (DSM-5) criteria. Disordered sleep was identified through parent survey and chart review by research staff. Participants with disordered sleep had either: 1) parent reported difficulty with sleep initiation or sleep maintenance; 2) ICD-10 diagnosis of disordered sleep (G47 or F51); or 3) a prescription for melatonin, clonidine, or mirtazapine with indication as a sleep aid. There was no difference in mean collection time between ASD subjects with (12:30PM ± 2:48) and without (1:00PM ± 3:00) disordered sleep (p = 0.34). The sleep disorder group was 18% female (14/77) and had a mean age of 56 (±16) months. The non-sleep disorder group was 14% female (9/63) and had a mean age of 56 (± 13) months.

Saliva collection and processing
Before collecting saliva samples, each subject rinsed their mouth with tap water. Approximately 1 mL of saliva was obtained through swab collection using an Oragene RNA collection kit (DNA Genotek; Ottawa, Canada). Samples were stored at room temperature until processing. A Trizol method was used to purify the salivary RNA and a second round of purification was followed using an RNEasy mini column (Qiagen). Yield and quality of the RNA samples was assessed with the Agilent Bioanalyzer. This was done prior to library construction in accordance to the Illumina TruSeq Small RNA Sample Prep protocol (Illumina; San Diego, California). Identification and quantification of saliva miRNA and microbial RNA was performed using next generation sequencing (NGS) on a NextSeq 500 instrument (Illumina), following the TruSeq Small RNA Library Preparation Kit protocol (Illumina, San Diego, CA). Alignment of mature miRNA reads was performed with the miRbase21 database using the Shrimp2 algorithm in Partek Flow software (Partek, Inc., St. Louis, MO). Mapping of unique microbial transcripts was performed using the K-Slam database, which references the NCBI Taxonomy database [46]. Taxons were defined by their family, genus, species, and subspecies (when available). The human miRNAs and microbial RNAs present in raw counts of 10 or more in at least 10% of samples were interrogated for oscillating expression. A quantile normalization technique was applied to the human miRNA and microbial RNA datasets separately, prior to statistical analysis.

Identification of oscillating salivary RNAs
A two-way analysis of variance (ANOVA) was performed using sample sets 1 and 2 based on binning the samples into their approximately replicated collection times, to identify host miR-NAs and microbial RNAs that varied significantly (FDR<0.05) with collection time but not the day of collection (in order to eliminate RNAs which could be influenced by daily variations in routine). A subset of miRNAs and microbial RNAs that were highly associated with time of collection (R ! 0.90 or 0.84 in sample sets 1 and 2, respectively; p<0.001) were then used in a naïve hold-out set (sample set 3) to assess predictive accuracy for time of collection with a multivariate regression analysis. The miRNAs that showed the strongest circadian oscillations were termed CircaMiRs and the microbes that displayed the strongest oscillations in transcriptional activity were termed CircaMicrobes. Relationships between CircaMiRs and CircaMicrobes were investigated with a Pearson Correlation analysis.

Functional Interrogation of CircaMiRs and CircaMicrobes
Classification of the mapped microbial RNAs within defined metabolic and functional categories was performed through conversion of microbial reads to Kyoto Encyclopedia of Genes and Genomes (KEGG) Orthology identifiers (KO IDs) which were mapped by MicrobiomeAnalyst software [47] to a set of 202 different KEGG Modules, KEGG Pathways, and COG Categories. For each participant, KEGG and COG data were summed across four collection periods (i.e. 7-9 AM, 10 AM-2 PM, 3-6 PM, 7-10 PM) for all of the days saliva samples were collected. Changes in expression of individual functional clusters were explored with a nonparametric Kruskal Wallis ANOVA. Patterns in functional clusters across the four time periods were visualized with hierarchical clustering analysis and a partial least squares discriminant analysis in MetaboAnalyst software.
The potential biologic impact of the CircaMiRs was investigated through functional annotation of their high confidence mRNA targets (p<0.05, Micro-T Score ! 0.95) in DIANA miR-Path v3 software and Ingenuity Pathway Analyst software (IPA, Qiagen). KEGG pathways over-represented by these mRNA targets were determined with Fisher's Exact test with FDR correction (FDR<0.05). Inter-relatedness of protein products for the mRNA targets was explored in String v10.5. Alignment of salivary RNA to the RefSeq Transcripts database in Partek Flow permitted quantification of local (oropharyngeal) mRNA targets for salivary Circa-MiRs (that were 50 base pairs). Relationships between CircaMiRs and mRNA targets were explored with Pearson's correlations.
To further explore the potential biological significance of the miRNA data, we examined the levels of the oscillating salivary CircaMiRs in the same cohort of 2-6 year old children with ASD examined for miRNA expression who either had normal sleep patterns (n = 63) or disordered sleep symptoms (n = 77). Group differences in mean salivary CircaMiR expression between the sleep disorder and non-sleep disorder groups were identified with a non-parametric Mann Whitney U-test. A two-way ANOVA assessed relationships between CircaMiRs, disordered sleep, and collection time, as well as sleep disorder-time interactions. Finally, a multivariate linear regression was used to determine the ability of the most robust CircaMiRs to predict collection time in the ASD children with and without sleep disorders.

Influence of daily routines on the oral transcriptome
To investigate the potential impact of daily routines on salivary miRNA and microbial RNA levels we examined associations between the oral transcriptome in sample set 3 and: 1) time since last meal (in hours); 2) time since last tooth brushing (in hours); and 3) time since last sleep (in hours). Significant relationships (|R|>0.40; FDR<0.05) between these three variables and salivary RNA levels were reported.
Supporting information S1