Longitudinal saliva omics responses to immune perturbation: a case study

Saliva omics has immense potential for non-invasive diagnostics, including monitoring very young or elderly populations, or individuals in remote locations. In this study, multiple saliva omics from an individual were monitored over three periods (100 timepoints) involving: (1) hourly sampling over 24 h without intervention, (2) hourly sampling over 24 h including immune system activation using the standard 23-valent pneumococcal polysaccharide vaccine, (3) daily sampling for 33 days profiling the post-vaccination response. At each timepoint total saliva transcriptome and proteome, and small RNA from salivary extracellular vesicles were profiled, including mRNA, miRNA, piRNA and bacterial RNA. The two 24-h periods were used in a paired analysis to remove daily variation and reveal vaccination responses. Over 18,000 omics longitudinal series had statistically significant temporal trends compared to a healthy baseline. Various immune response and regulation pathways were activated following vaccination, including interferon and cytokine signaling, and MHC antigen presentation. Immune response timeframes were concordant with innate and adaptive immunity development, and coincided with vaccination and reported fever. Overall, mRNA results appeared more specific and sensitive (timewise) to vaccination compared to other omics. The results suggest saliva omics can be consistently assessed for non-invasive personalized monitoring and immune response diagnostics.

Precision medicine continues its rapid development toward clinical applications aided by new sequencing technology and computational capabilities. Major efforts have concentrated on evaluating disease risk from genomic information 1,2 , including direct to consumer platforms, like 23andMe 3 , as well as pharmacogenomic evaluations 4 . Implementing omics profiling in the clinic will require evaluation of patients over time. The utility of such profiling has been evaluated in individual monitoring pioneered in the integrative Personal Omics Profiling (iPOP) study 5 , and expanded recently to include profiling using electronic health devices 6 . The Pioneer study 7 also incorporated behavioral coaching to improve clinical biomarkers based on participants' individual data. Additional developments have included utilizing host-microbiome data in insulin resistant individuals in a study of weight gain 8 and in prediabetics 9 , investigating biological age 10,11 as well as monitoring of astronauts in the recent NASA twin study 12 .
In this investigation we are extending integrative omics to evaluate the utility of such monitoring using saliva. There has been long-standing interest in saliva for non-invasive diagnostics and health monitoring, and saliva omics is an emerging field, with broad profiling that includes total saliva RNA and proteomes, as well as cell-free RNA identification, extracellular vesicle (EV) profiling, miRNAs as biomarkers, and salivary microbiomes [13][14][15][16][17][18][19][20][21][22][23][24][25] . Utilizing saliva for non-invasive monitoring is important in evaluating vulnerable populations, including infants, children, older adults and immunocompromised individuals. Additionally, saliva is important in evaluation of health in remote or underserved locations, when limited resources are available, where processing of blood samples might not be feasible, or a physician may not even be available. Such monitoring is also of particular interest for evaluating active personnel, including astronauts in deep space missions. The recent twin astronaut study evaluated multi-omics utility but also highlighted the logistic issues of using blood samples when these cannot be processed on-site 12 . The COVID-19 pandemic has additionally ignited interest in the use of saliva for rapid diagnostics, towards a rapid and minimally invasive diagnostic that can be used without risk to personnel www.nature.com/scientificreports/ (possibly as a home-use kit), including profiling viral loads (from posterior oropharyngeal samples) 26 , and current work continues to evaluate the sensitivity of saliva for practical implementations [27][28][29][30][31] .
We are carrying out a clinical trial monitoring individualized response to pneumococcal vaccination, and in a proof-of-principle case-study, we monitored individualized response to the standard 23-valent pneumococcal polysaccharide vaccination (PPSV23), in a generally healthy individual (Caucasian male, 38, has reported chronic sinusitis), and carried out integrative profiling on saliva pre-and post-vaccination with pneumococcal PPSV23 vaccine. This is to our knowledge the most extensive saliva-focused omics dataset on an individual, covering 104 timepoints over one year. The period covers a healthy period as well as monitoring of innate and adaptive immune responses following vaccination. Protein and RNA from saliva were produced over 100 timepoints over the course of 1 year, and comprehensive untargeted proteomics and RNA-sequencing were carried out for all samples. The saliva sampled timepoints included three periods of particular reference in this manuscript: (1) 24 h hourly sampling without intervention to assess a healthy hourly baseline, (2) 24 h hourly sampling that included vaccination with pneumococcal vaccine (PPSV23) to assess response to the vaccine, (3) daily sampling following the vaccination to assess potential innate and adaptive immune responses reflected in the molecular saliva components.
Our study reveals multiple changes in response to pneumococcal vaccination that are observable in saliva. The microscopic collective behavior of multiple omics reflects physiological changes associated with immune response, including fever, innate and adaptive responses profiled over multiple scales. This case study provides a resource for future saliva studies, towards more effective non-invasive diagnostics.

Results
Samples and assays. We followed a single individual (m, 38, Caucasian), in general good health (has reported chronic sinusitis) over the span of a year. To observe whether the effects of perturbation can be profiled in saliva we carried out the profiling over 3 time frames. In the first 24 h time frame (TFH1) we established a baseline, obtaining a saliva sample from the subject hourly without perturbation over his standard routine. In the second 24 h time frame (TFH2), the subject was vaccinated with pneumococcal polysaccharide vaccine (PPSV23) within 3.5 h of waking up (at 10.30 am), while otherwise maintaining a similar routine as in the first period (including food intake and meal timing), and again saliva samples were taken hourly. We should note here that the subject reported fever ∼ 7.5 h post the vaccination (between timepoints at 5 and 6 pm), lasting for about 4 h (10 pm). The two time periods, TFH1 and TFH2, were treated as paired and combined in the analysis below ( TF ) to identify changes induced by the vaccination, by effectively removing daily normal routine effects for this individual. Additionally, in the third time frame (TFD) we monitored the subject daily for over a month, pre-and post-vaccination to identify potential immune changes over both innate and adaptive time frames (Fig. 1).
The daily samples were all taken at 8 am, to limit variability. Unstimulated saliva was collected both for downstream total RNA profiling, mass spectrometry proteomics, as well as for extraction of extracellular vesicles which were profiled for various small RNA molecular species (both host and non-host)-see "Methods" for further sample details. Figure 1. Study overview. The study followed an individual over the span of a year, collecting 100+ saliva samples. In this manuscript we discuss three time frames of interest for saliva sampling: (1) TFH1, where 24 hourly samples were taken, (2) TFH2, where 24 hourly samples were taken but the subject was also vaccinated with PPSV23 during this time frame, (3) TFD, where daily samples were taken. At each timepoint unstimulated saliva was collected (at a rate ∼ 500 µ l/min ). 3 ml of saliva were stored directly and subsequently used to extract extracellular vesicle RNA, and proteins for mass spectrometry proteome profiling. Furthermore, 2 ml saliva was collected with Oragene (DNAGenotek) kits, which contain a stabilizer, and used to profile total saliva transcriptomics using RNA-sequencing; see "Methods" for further collection and processing details. The data were used to generated time series with MathIOmica, which revealed multiple trends corresponding to response to immunization with PPSV23. Total saliva proteomics. We profiled the total saliva proteome, using isobaric tandem mass tags (TMT) for quantitation using LC-MS/MS (liquid chromatography followed by mass spectrometry). We identified 12,473 proteins overall, with 4141 proteins (UniProt identifiers 38 ) based on 2 unique peptides per protein, (11,005 proteins based on 1 unique peptide per protein) overall across all 95 samples where proteomics was carried out. Relative protein intensities were computed against a common pooled sample comprising of multiple healthy (pre-vaccination) weekly samples that was used across all TMT sample pools. The data were thus combined, and normalized using a Box-Cox 39 transformation to obtain normal distributions. To construct the time series, the data were filtered again for 2 unique peptides, having less than 1/4 missing values, and no constant time series to obtain 724, 956, 759, and 662 proteomics time series for TFH1, TFH2, TF and TFD respectively. All timepoint intensities were defined with respect to the first timepoint intensity for the hourly series for each respective protein, and to the vaccination day for TFD.

Analysis of saliva extracellular vesicles.
In addition to considering total saliva, we also implemented consistent extraction of EVs from 1 ml saliva, using ExoQuick-TC (SBI) and an overnight precipitation to obtain EV pellets, from which we extracted RNA. We carried out nanoparticle tracking analysis (ZetaView, Particle Metrix) and recorded median concentrations of 6.2 × 10 10 particles/ml with EV peak of 114.5 ± 4 nm (Fig. 2). The extracted EV RNA was sequenced (small RNA-sequencing) and the results were mapped to multiple databases using the exceRpt 40 through http://genbo ree.org 41 . These included GENCODE transcripts, PIWI interacting RNA (piRNA), micro RNA (miRNA), and exogenous genomes and exogenous ribosomal RNA (rRNA) genomes, and multiple other biotypes (Fig. 3). The various biotypes detected have different variabilities (Fig. 3a). The majority of detected reads are from exogenous genomes ( ∼ 10 6 ) and protein coding (GENCODE transcripts, ∼ 10 6 as well). The different biotypes per sample (Fig. 3b), for the TFH1, TFH2 and TFD time frames are indicative of a change in the relative distributions of the biotypes for the daily samples following the vaccination (increase of exogenous genomes content). This may be partly attributed to sampling as all TFD samples were taken at 8 am, and the corresponding early samples for TFH1 and TFH2 are also similar in biotype relative abundances, but differ in later samples during each day.
In terms of the exogenous genomes, taxonomy trees were constructed per sample, and also for the aggregate samples using Genboree 42 (Fig. 4a). The majority of abundances were assigned to bacteria (89.5%), and to Eukaryota (6.4%), where in terms of majority assignments at the next level, Bacteroidetes/Chlorobi group (28.5%),  www.nature.com/scientificreports/ 27.2% were assigned to Proteobacteria and 17.1% to Firmicutes. Clustering of the overall top taxa by normalized read count was indicative of consistency across samples, with no significant sub-grouping at this level (Fig. 4b).
We again carried out downstream analysis in MathIOmica 37 , to create EV time series for mapped data for each of the three time frames (TFH1, TFH2, TFD), and the paired difference ( TF ) for GENCODE mapped entities, host piRNA, host miRNA and exogenous rRNA and exogenous taxa. Time series were constructed for entities where 0 count values were tagged as Missing, counts < 1 were considered equivalent noise, transcripts with more than 1/4 timepoints missing were removed, and transcripts with constant values across all timepoints were also removed.
Temporal trends identified in saliva. Time series classification. The time series for all omics discussed below were classified into temporal trends using MathIOmica's 37 spectral methods that classify signals based on their autocorrelations, i.e. correlating a time signal with a delayed version of itself, where the delay is characterized as a time lag (e.g. lag 1 corresponds to a delay of 1 time interval unit). The method uses a Lomb- Scargle   a   b   IG_J_pseudogene  TR_D_gene  TR_J_pseudogene  IG_D_gene  vaultRNA  IG_J_gene  translated_unprocessed_pseudogene  sRNA  transcribed_unitary_pseudogene  IG_C_pseudogene  bidirectional_promoter_lncrna  TR_C_gene  ribozyme  non_coding  TR_J_gene  scaRNA  3prime_overlapping_ncrna  TR_V_pseudogene  IG_V_gene  IG_V_pseudogene  IG_C_gene  circularRNA  pseudogene  polymorphic_pseudogene  Mt_tRNA  non_stop_decay  Mt_rRNA  rRNA  unitary_pseudogene  snoRNA  TR_V_gene  macro_lncRNA  transcribed_processed_pseudogene  transcribed_unprocessed_pseudogene  sense_overlapping  sense_intronic  unprocessed_pseudogene  TEC  piRNA  exogenous_miRNA  processed_pseudogene  snRNA  antisense  lincRNA  nonsense_mediated_decay  miRNA  tRNA  exogenous_rRNA  processed_transcript  retained_intron  misc_RNA  protein_coding   www.nature.com/scientificreports/ transformation [43][44][45] to generate periodograms whose inverse Fourier transform can then produce a set of autocorrelations at different lags for a given time series. Our classification successively identifies time series from the dataset that have statistically significant autocorrelations at particular time lags. In summary (see "Methods"), MathIOmica's classification generates three sets of classes, strictly based on temporal behavior: (1) significant autocorrelations at various lags, (2) no autocorrelations, but with positive spikes (abnormally high signals above baselines present at single timepoints), (3) no autocorrelations, but negative spikes present (abnormally low signals below baseline at single timepoints). Within each class a two-tier classification into groups and subgroups  www.nature.com/scientificreports/ is carried out: this approach first separates within-class autocorrelation groups by clustering on autocorrelation lags: signals that may have statistically significant autocorrelation for the class lag, but may still exhibit underlying different structure at other lags. Additionally, the second level clustering into subgroups is based on intensities, and allows us to separate signals that may have different phase (directionality/sense), which cannot be obtained from the periodograms. The analysis was carried out for each of the omics individually and thousands of individual component trends were identified in the different classes and subgroupings therein. A brief summary is provided in Table 1. The entirety of visualizations and classification memberships are available in the online data files (ODFs), including heatmaps per omic per individual time frame trends, as well as all the code to generate these. We also combined all the classified information to obtain an integrated view of the various omics. Below we showcase parts of the mRNA analysis, as well the results of all the omics combined.
Saliva mRNA data analysis. The trends shown in Fig. 5 correspond to the mRNA time series showing statistically significant time series trends (p value < 0.01 based on bootstrap simulations, n = 100,000) for each of the time frames, for Lag 1 classes.
Hourly results (TFH1, TFH2 and TF ). The saliva mRNA showed variation across the day in the untreated TFH1 period. Overall 5781 time series of mRNA isoforms were found to have statistically significant trends, with 1085 Lag 1, 6 Spike Max, 3597 Spike Min and 1093 other Lag class memberships. The Lag 1 group is shown in Fig. 5a, where the 1085 time series are further assigned into groups and subgroups based on clustering, according to their different temporal behaviors as described above. In the Lag 1 classification in Fig. 5a there are two groups (G1 and G2). G1 has 3 subgroups (S), S1, S2 and S3 with 187, 522 and 373 time series in each respectively. G2 has 2 subgroups S1, S2 with two and one time series respectively. The groupings show substantial variation in these isoforms' intensities during the day, with G1S1 showing gradual decreases, G1S2 peaking after morning until the evening, and G1S3 showing peaks later in the evening and night (the first timepoint is at 7 am).
In the analysis of the 24 h period spanning vaccination, TFH2, we should note that the subject reported fever  www.nature.com/scientificreports/ Given the variation observed in TFH1, we constructed the TFH time series, using paired differences of intensities at each timepoint. The approach aimed to remove non-vaccination daily variation, and resulted in 1200 time series with statistically significant trends. These included 369 Lag 1, 3 Spike Max, 2 Spike Min, and 826 other Lags memberships. The Lag 1 results are shown in Fig. 5c. The subgroupings of 273 G1S1 isoform time series, show punctuated trends following vaccination and also coincidental with the reported fever, lasting about 4 h (timepoints). Furthermore, the G1S2 subgroup of 94 time series is indicative of up-regulation following the vaccination. Additionally a distinct up-regulation of a subset of genes is observed to coincide with the reported fever (Fig. 5c).
We carried out Gene Ontology (GO) 47 and Reactome Pathway enrichment analysis 46,48 and identified multiple involved pathways. Results for TF with False Discovery Rate, FDR < 3 × 10 −3 are shown in Table 2, and full results available in the ODFs. For the set of TF genes showing immediate response post vaccination and response during the fever period (Class Lag 1, G1S1 Fig. 5c), results include (Table 2A (i)) endosomal/vacuolar pathway, antigen presentation (class 1 MHC) and processing, interferon gamma and alpha/beta signaling, neutrophil degranulation and ER-Phagosome pathways indicative of the immune activation. Furthermore, a set of genes that show continued up-regulation following the vaccination (Class Lag 1, G1S2 Fig. 5c) had enrichment of various immunological pathways including TCR signaling-related pathways, PD-1 signaling, also Interferon gamma signaling, Costimulation by the CD28 family, MHC class II antigen presentation, Cytokine Signaling in Immune system, and also Neutrophil degranulation pathways and general Immune System pathways.
Daily results (TFD). We also monitored the daily changes post vaccination for 1 month. For the mRNA analysis, 3762 time series had statistically significant trends. These included 439 Lag 1, 2 Spike Max, 2248 Spike Min, and 1073 other Lag memberships. The Lag 1 TFD results are shown in Fig. 5d. As shown in the figure, in the G1S1 subgroup 219 isoform time series show up-regulation, for about 11 days post vaccination, followed by a return to lower expression levels. The 218 time series in the G1S2 subgroup also show a later up-regulation response, after 11 days compared to G1S1, again following the vaccination, and lasting for the remainder of the daily observation period.
Reactome pathway and GO enrichment analysis also identified multiple pathways corresponding to each trend. Reactome results for TFD Lag 1 are shown in Table 2B(i)-(ii) for Lag 1 G1S1 and G1S2 subgroups (full results, including GO terms available in the ODFs). In the TFD Lag 1 G1S1 group, over-representation included Endosomal/Vacuolar pathway (16 genes) Interferon Signaling (29 genes), Cytokine Signaling in Immune system (59 genes), Antigen Presentation (MHC I related, 14 genes) and inteleukin 4 and 13 signaling pathways (18 genes), and Immune System (97 genes). These are indicative of an early response within the first days after The TF results correspond to paired differences between TFH2 and TF1 hourly points, to remove intra-day variation so as to focus on the perturbation vaccination responses. The plot is indicative of a response to the vaccine and a response that coincides with the reported fever. (d) For the daily data, TFD, the corresponding vaccine timepoint is Day 3. There is a direct response the days following the vaccination, and a different response in a subset of genes approximately a week following the vaccination, corresponding to immune system activation. www.nature.com/scientificreports/ Cytokine signaling in immune system 33 Insulin-like growth factor-2 mRNA Binding Proteins (IGF2BPs/IMPs/VICKZs) bind RNA 5 1.6 × 10 −7 1.0 × 10 −5 www.nature.com/scientificreports/ vaccination. For TFD Lag 1 G1S2 the results included general Immune System activation (88 genes), and also Cytokine Signaling in Immune System (51 genes) as pathways with statically significant over-representation and having the most genes identified.
Other omics. In addition to the mRNA, each other set of omics was individualized analyzed to identify temporal trends, using the same classification method in MathIOmica as described above. This identified statistically significant trends for time series for different omics in the different time frames are shown in Table 1. The different classes for the omics datasets were joined within each respective time frame, and data within each combined class were clustered together. The breakdown of identified trends included overall 35,372 time series for TFH1, 25,739 for TFH2, 18,815 for TF , and 27,094 for TFD. In reference to the corresponding omics, the EV GENCODE identifiers accounted for more that 78% of the time series across all time frames. The results from Lag 1 are again shown in Fig. 6. In terms of temporal behavior, we notice again similar responses in the various time frames. In TFH1, we again see large temporal variation, with sets of time series being up-regulated during awake time and a subset during night time as shown in Fig. 6a. In TFH2 (Fig. 6b), the effects of the vaccination become apparent with up-regulated responses following the vaccination.
In the paired comparison TF for Lag 1 in Fig. 6c, there is a major subgroup (G1S1) with 7700 time series. A subset of the omics show increases in intensity about 2-3 h post vaccination for a few hours, and additionally again increase in intensity about 14 h post vaccination (after the fever time span). We note here that the response/time pattern appear lagging compared to the corresponding mRNA results in Fig. 5c, by approximately 3 h, both following the immediate vaccine response, and also for the reported fever time span. The mRNA over-representation analysis was discussed above. Proteomics displayed similar trends to the mRNA. For Lag 1 the Reactome analysis for the aggregate group/subgroup proteins resulted in multiple statistically significant pathways ( FDR < 0.01 ), where the top pathways (FDR < 5.9 × 10 −13 ) included Neutrophil degranulation (96 entities identified), Innate Immune System (135 entities), Immune System (with 158 entities identified), and more (the top 20 ranked by smallest FDR are shown in Table 3A, FDR < 2.6 × 10 −7 ).
GO analysis for the EV GENCODE results for TF Lag 1 were generic and included multiple cellular processes. The results appear as non-specific to immune response, particularly as no over-representation in Reactome pathways was found to be statistically significant based on FDR. For EV miRNA the TF Lag 1 showed statistically significant over-representation in GTP binding (p value < 4.6 × 10 −14 ), perinuclear region of cytoplasm (p value < 9.3 × 10 −14 ), nerve growth factor receptor signaling pathway (p alue < 1.5 × 10 −13 ), MAPK cascade (p The TF results correspond to paired differences between TFH2 and TFH1 hourly points, to remove intra-day variation so as to focus on the perturbation vaccination responses. The plot is indicative of a phased response to the vaccine compared to the mRNA responses and a response that is again shifted compared to the reported fever (cf. Fig. 5c). (d) For the daily data, TFD, the corresponding vaccine timepoint is Day 3. There is a direct response the days following the vaccination. The EV omics dominate the information in this plot, compared to the mRNA in total saliva response (cf. Fig. 5d).   www.nature.com/scientificreports/ at lower levels, and a set that remains constant in intensity for the duration of the month's measurements and increases towards the end of the monthly period. The mRNA TFD pathway analysis results were discussed above (Table 2B). Proteomics time series on the other hand displayed fewer pathway enrichment results as shown in Table 3B, including various Nonsense-Mediated Decay pathways, Translation related pathways, ROBO receptor signaling and others. GO analysis for the EV GENCODE results for TFD Lag 1 included generic molecular functions such as protein binding and ATP binding and biological processes relating to protein phosphorylation (105 IDs, adjusted p value < 4.7 × 10 −5 ). Full enrichment analysis results for all classifications and omics are available in the ODFs on Zenodo.

Discussion
We have presented here our findings from a case study of the utility of saliva towards personalized health monitoring. Following vaccination of a subject with pneumococcal vaccine (PPSV23) we were able to detect distinct signatures in various saliva omics. We were able to profile more than 65,000 components in various time frames over time, and identify 18,000+ time series that had statistically significant temporal trends. The time series trends observed were indicative of immune response, which coincided in timing with the vaccination, and fever reported by our subject. The time frames of immune responses observed are concordant with our expectations of innate and adaptive immunity development, as seen both in the immediate hourly as well as the short-and long-term daily responses observed. Various pathways were activated, involved in immune response and regulation, including interferon signaling and MHC antigen presentation. The immune activation spanned an initial response within hours, as well as long term response extending for over a month.
Our results suggest that saliva omics can be consistently assessed for personalized monitoring. While multiple omics provide responses post-vaccination as discussed for the TF and TFD time frames, mRNA results appear more specific and sensitive (timewise) to the vaccination. EV response, particularly transcript-level (GENCODE), though very sensitive and responsive does not yield as specific results as the mRNA from total saliva in a nontargeted approach. In fact many EV mRNAs show statistically significant temporal responses, which is likely due to increased EV release from various cells, but not specific in terms of reflecting the functional responses in the cells of origin. EV results also suggest a lagged response by a few hours compared to the mRNA observations, which also suggests that the mRNA measurements can potentially provide more timely data for practical health monitoring. Additional omics showed responses concordant with the mRNA responses, including miRNA, piRNA and exogenous taxa quantitations, however these need further validation, particularly as our knowledgebase of pathway implications and functional important of such omics is still under development. Some of the omics time series are found in responses across hourly and daily samplings (Table 4), and such sets of omics can be targeted for non-invasive health monitoring. The processing of samples for mRNA involved the use of standardized kits that can be used by subjects remotely, and can facilitate storage of samples for about a month without refrigeration. Additionally, the mRNA sequencing preparation and result analysis are considerably faster than the other omics processing, so our recommendation based on our findings is to utilize similar approaches for mRNA broad profiling, while using targeted protein/EV-content assays. These omics must be coupled with standard physiological and molecular measures already utilized in the clinic for a complete assessment of health status. While our results are based on PPSV23 activation, we anticipate that they can be extended to additional vaccine and immune disease profiling, particularly with the goal of discovering immune-specific signatures for each affliction and/or intervention.
Previous work on saliva omics has focused on the evaluation of omics as biomarkers, with less emphasis on temporal changes, including proteomic and transcriptomic evaluations, EV characterization, miRNA profiling and microbiome mapping [13][14][15][16][17][18][19][20][21][23][24][25]49,50 53 . Such efforts were substantially improved by Grassl et al. who used deep profiling by mass spectrometry to identify more than 3700 human proteins 54 . Their study also assessed intra-day changes for waking and postprandial collection times, identifying enrichment of proteins associated with 'antibiotic' keywords in waking versus postprandial collection times. Our study significantly extends the temporal profiling aspect to hourly 24 h period monitoring (for 4141 proteins-UniProt IDs, with > 2 unique peptides used per identification), in addition to daily data for a month. We also observed substantial variability during the course of a day. To account for this daily variation, we used a paired two 24 h period comparison to elucidate changes particular to the vaccination (essentially subtracting normal hourly variation effects), and also limited collection to a single morning timepoint in our daily collection. Other studies of the saliva proteome have included the integration of transcriptomics with antibody-based proteomics [55][56][57] to assess the salivary gland content (Human Protein Atlas), suggesting up to 15,218 proteins are expressed in the gland itself, though these studies did not dynamically profile secreted saliva.
Salivary proteomics (at various scales) has been used in prospective clinical applications ranging from periodontitis, oral and other cancers, diabetes, Sjogren's syndrome, and to assess viral proteins in Zika virus, Dengue virus, HBV and HCV (review by Katsani and Sakellari 58 ). Saliva immunoglobulins levels in COVID-19 were also evaluated by Isho et al. 59 . With respect to pneumonia, Klein Kremer et al. measured overall increases in aggregate salivary protein levels in children diagnosed with Lobar pneumonia 60 . Recently Tsai et al. reported immunoassay results on 9 cytokines and C-reactive protein (CRP), and detected higher levels of CRP and IL-6 in children with pneumonia 61 . In our MS-based study we did not detect CRP/IL-6 changes directly, but we identified multiple proteins associated with immune pathways (Table 3) www.nature.com/scientificreports/ markers has been reviewed 62 , MS-based proteomics evaluation of longitudinal saliva responses in pneumonia and PPSV23 vaccination have not been carried out prior to our study, to the best of our knowledge. In addition to the coupled transcriptome-proteome evaluations (Human Protein Atlas, see above), focused saliva transcriptomics have also been previously evaluated, including through expression microarray analysis 15 , and high throughput sequencing by Spielman et al. who detected the expression of > 4000 coding and noncoding genes 17 . In our RNA-seq results we detected 67,319 GENCODE annotation transcripts showing non-zero values for at least 3 timepoints (81,098 observed in at least 1 timepoint), with more than 7493 transcripts detected consistently over 3/4 of the hourly observations, and 8155 over 3/4 of the daily observations.
We have also profiled the RNA of salivary EVs. There has been expansive interest in the diverse RNA content of EVs 63 , including by the Extracellular RNA Communication (ERC) program 64 , for applications in liquid biopsies, as markers in disease states and for cell-free precision medicine diagnostics. EVs are being evaluated as mediators of intercellular communication through molecular transport, offer stable containment of RNA, and can easily be collected for potential diagnostics 63 . EVs have been detected and may move across biofluids, with RNAs from bacteria, fungi, and other species having been reported in human plasma and saliva [65][66][67][68] 68 , processing an average 16 ×10 6 reads and mapping across various biotypes [with ∼ 34% mapping to hg19, 0.02% to piRNA, 1.65% to mature miRNA, including the identification of 336 miRNAs (detected in at least one out of 46 samples, with > 10 counts; 149 miRNAs in at least 50% of the study samples), and a large amount of reads not mapping to human transcriptome]. Godoy et al. compared EV RNA contents across multiple biofluids, including parotid saliva and submandibular and sublingual saliva 67 , and also detected low mapping to human transcriptome (as expected given the microbiome content in saliva), and observed 395+ miRNA ( ≥ 10 reads per million)and < 0.01% piRNA in their saliva samples, while also detecting multiple maps to GENCODE annotations, and other small RNA subtypes. In our mapping, we adopted the same mapping strategy as the exRNA Atlas, implemented by Rozowsky et al. 40 , with similar multi-biotype mapping results, including to non-human exogenous taxa, and our time series included 140-258 miRNAs, 275-589 piRNAs, and 55,499-58,863 GENCODE transcripts that could consistently be evaluated over time (3/4 of the samples in each time frame).
While the content of saliva EVs has been explored, their longitudinal changes, and in particular in response to pneumococcal vaccination (or pneucoccal disease) had not been investigated prior to our study. Additionally, the microbiome EV content is an area of new study, that has not been fully evaluated for its effects in pneumonococcal disease, and there is considerable interest in bacterial EVs (BEVs), for example in the context of cancer 73 . Our goal in this study was to focus on the host, so we did not evaluate the oral microbiome, beyond EV content, though this has been extensively studied, particularly with 16S ribosomal RNA profiling 74,75 , as previously reviewed 76,77 . Recent studies have included longitudinal monitoring of the oral microbiome in the context of oral health: Dzidic et al. investigated the long term effects of colonization during development as associated with tooth decay, associated carries with temporally divergent microbial constitution 78 . Kennedy et al. investigated oral microbiota composition using sequencing in children sampled at 6, 12 and 24 months of age 79 . Kahharov et al. longitudinally profiled the oral microbiome maturation of the Oral Microbiome in caries-free children 80 . Lif Holgerson et al. recently studied the longitudinal development of salivary microbiome in adolescents 81 . In future investigations it will be interesting to study further the interplay between host transcriptome/proteome and oral microbiome in infectious disease, and monitor these in parallel, as we do expect the microbiome to directly affect EV content, and partake in multi-omic interactions potentially modulating immune responses.
In terms of the longitudinal monitoring of individuals over time, profiling multiple omics, such approaches were pioneered with the iPOP study 5 , that measured up to 20 timepoints of multiple blood-based omics in an individual, over a timeframe that included coincidental viral infections, and the onset of type 2 diabetes. David et al. 82 monitored the daily microbiome in gut and saliva of two individuals, with 274 saliva samples profiled for 16S ribosomal RNA, with their findings indicating that travel and enteric infections affecting community structure. The Pioneer study 7 which incorporated behavioral coaching to improve clinical biomarkers based on 108 participants' individual data, including bood-based multiomics profiling, also measured saliva cortisol and dehydroepiandrosterone (dhea) levels for stress assessment every three months. Additional longitudinal studies have focused on host-microbiome characterization of multiple insulin resistant individuals and studying weight gain 8 and in prediabetics 9 , investigating biological age 10,11 as well as monitoring of astronauts in the recent NASA twin study 12 . Our study extends these approaches, not only by longitudinally monitoring saliva host transcriptome, proteome and EV content simultaneously, but also in providing dense profiling with hourly sampling. Our study, in conjunction with the previous individual monitoring efforts provide the first steps in personalized wellness monitoring, not only in demonstrating the feasibility of utilizing state-of-the art multiomics technologies, but also providing extensive datasets for modeling temporal processes in direct applications to health, including in our case non-invasive monitoring of immune responses, such as vaccination.
Our study also has limitations: even though we attempted to pair time responses for the hourly data, this is still a single subject case-study (n = 1), and our results will need to be validated. Furthermore, due to limited samples and resources, we did not carry out immune profiling, such as cytokine assays or functional assays to assess the immune acquisition to the different components in the PPSV23 vaccine. We were also unable to obtain blood samples across all the time points, as our focus was a first evaluation of saliva omics. Additionally, our study did not specifically target the salivary microbiome, and the meta-analysis of EV RNA content indicated substantial www.nature.com/scientificreports/ variability in overall taxa, the composition of which is expected to vary across individuals. In the next stage of our long-term project we have already collected samples from multiple subjects being vaccinated at the same time points with PPSV23 and monitored over multiple timepoints. Given further resources, our goal is to utilize these samples to both validate and extend our findings to also include monitoring of blood components. By comparing the responses in blood and saliva we will be able to assess to what extend saliva may be used as a proxy for blood monitoring, identifying common and different responses in different tissues. Finally, in monitoring total saliva omics in bulk, we are ignoring the multi-cellular composition of saliva. With the availability of single-cell RNA-seq methodology, we anticipate that we will be able to also assess the cell-type-specific response in saliva. In summary, saliva provides a promising venue for non-invasive diagnostics of immune response. This is particularly important for enhancing our diagnostic capabilities for multiple viral or bacterial responses, particularly in cases where blood may not be easily available, due to technical issues, remote locations (e.g. monitoring active personnel), lack of specialized equipment and healthcare availability (e.g. due to socio-economic factors), patient vulnerability (immuno-compromised, children, and elderly populations). Given the current pandemic (COVID-19), enhancing our diagnostic capabilities has become a high priority. While the utility of saliva for differentiating between different afflictions still needs to be evaluated, our study provides the first steps towards a no-pain no-blood diagnostic process that can greatly enhance our capabilities for universal individualized health care and diagnostics.

Methods
Data and protocol availability. Sequencing data reported here are available on Gene Expression Omnibus under accessions GSE108664 (Saliva mRNA-sequencing) and GSE108666 (EV small RNA). Proteomics data are available on MassIVE as part of accession MSV000081869. All scripts and data analysis code utilized in the integrative analysis are available on Zenodo (DOI:10.5281/zenodo.3987587) as online data files (ODFs), in addition to results and methods as referred to in the manuscript.
Sample collection. Samples were taken in the three time frames hourly for TFH1 and TFH2, and daily for TFD. In TFH2, the PPSV23 pneumococcal vaccine was administered at approximately 10.30 am (between sample collections at 10 and 11 am). Following the vaccination, and after the 24 h monitoring, daily samples were taken for about a month. At each timepoint 5 ml saliva were collected always in the same order: 2 ml in an Oragene (DNAGenotek) tube for RNA sequencing, and 3 ml in a conical tube for EV characterization and mass spectrometry proteomics, as described below. The collected samples were from unstimulated saliva (passive drooling), where the subject was instructed to let saliva collect in the front of their mouth and spit as saliva accumulated over time. The collection took approximately 10 min total per timepoint, for an estimated 500 µ l/min unstimulated salivary flow rate for the subject. Conical tube samples were immediately stored in a −20 • C noncommercial freezer, prior to transfer to the laboratory on ice where they were immediately stored at −80 • C on receipt. Oragene tube samples were capped and mixed with the stabilizing liquid that is part of the Oragene tube, and then kept at room temperature until transfer to the lab, where they were processed for RNA-sequencing as described below. Daily samples were all taken at 8 am, to limit variability. Additionally the subject followed the exact same diet and meal timings during TFH1 and TFH2, and had neither meals/drinks nor teeth brushing for at least 30 min prior to TFH1 and TFH2, and 1 h prior to TFD sample donations.
Saliva sample processing for RNA-sequencing. The saliva samples (2 ml) for RNA processing were collected in Oragene (DNAGenotek) tubes. The samples were incubated at 50 • C for 1 h, and stored at −80 • C . RNA Processing: 500 µl aliquots were incubated at 90 • C for 15 min in a heating block and then then cooled to room temperature. 20 µl neutralizer solution were mixed with each aliquot, vortexed and incubated on ice for 10 min, precipitating impurities and inhibitors. The sample was then centrifuged at 13,000 g for 3 min. The supernatant was transferred into a new microcentrifuge tube, 2 volumes of cold 95% EtOH were added and mixed thoroughly, followed by incubation at − 20 • C for 30 min. Following centrifugation at 13,000g for 3 min, the precipitate was collected. This pellet was dissolved in 350 µl of buffer RLT (RNeasy Micro kit). 350 µl of 70% ethanol were added and mixed. Additional steps followed the RNeasy Cleanup (Qiagen) per manufacturer's instructions to obtain concentrated RNA.
Libraries were constructed and sequenced by Novogene, using a Eukaryotic directional mRNA library (NEB). cDNA preliminary concentration was quantitated on a Qubit (Life Technologies), an Agilent 2100 Bioanalyzer was used to test the insert size, and Q-PCR was used to quantify the library effective concentration precisely. The cDNA libraries were sequenced on an Illumina HiSeq 4000.
Saliva EV processing. Saliva samples for EV processing were collected in a conical tube and stored in − 80 • C on sample receipt. EVs were processed from 500 µl saliva, following centrifugation at 3000g for 20 min at 4 • C . 500 µl of saliva were centrifuged at 3000g for 20 min at 4 • C to remove cells and cell debris. ExoQuick-TC Exosome Precipitation Solution (SBI) was added to the supernatant in a 2:1 ratio, and the mixture was refrigerated overnight at 4 • C . Following incubation samples were centrifuged 1500g for 30 min at 4 • C . The supernatant was aspirated and residual ExoQuick-TC solution was spun down at 1500g for 5 min. EV pellets were stored at −80 • C . RNA was extracted using the SeraMir RNA Amplification Kit (SBI) per manufacturer's instructions. The quality of EV RNA was checked using a 2100 Bioanalyzer (Agilent).
Small RNA library preparation was carried out using NEBNext Multiplex Small RNA library prep kit (New England Biolabs) following manufacturer's instructions. After PCR amplification, quality of libraries was assessed using a high sensitivity DNA kit on a Bioanalyzer (Agilent) according to manufacturer's instructions. Size selection was performed using 3% agarose dye-free marker H cassettes on a Pippin Prep (Sage Science) following www.nature.com/scientificreports/ manufacturer's instructions with a specified collection size range of 125-153 bp. Libraries were further purified and concentrated by ethanol precipitation, resuspended in 10 µl of 10 mM tris-HCl (pH = 8.5) and quantified using a Qubit and a Bioanalyzer. Based on the quantification, equimolar library pools were prepared, quality was assessed as described above and the library was further diluted to 4 nM using 10 mM tris-HCl (pH = 8.5). Pooled libraries were sequenced at a final concentration of 1.2 pM on an Illumina HiSeq 2500 (15-plex, 1 × 50 bp format).

Exosome quantitation by ELISA. EV concentrations were quantitated using the EXOCET Exosome
Quantitation Assay kit (SBI). EVs from 1 ml of saliva were precipitated using the Exoquick TC protocol (see above). Each exosome pellet was dissolved in 80 µl of lysis buffer and diluted with 80 µl of PBS to be used for duplicate reactions. Samples were then incubated at 37 • C for 5 min to liberate EV proteins, vortexed for 15 s, and centrifuged at 1500g for 5 min to remove debris. The supernatant EV protein samples were then assayed on a microtiter plate following the EXOCET kit manufacturer protocol (SBI), including 7 standards and blanks in duplicates. The plate was read using a spectrophotometric plate reader (Bio-RAD) at 405 nm. Spectrophotometry results for standards were used to obtain a linear fit, and sample results were indicative of ∼ 10 9 EVs/ml (supplementary data on Zenodo).
EV transmission electron microscopy. Isolated EVs were fixed in 2% paraformaldehyde (PFA) for 5 min. For negative-staining of EVs, 5 µl of the sample solution was placed on a carbon-coated EM grid and EVs were immobilized for 1 min. Next, the grid was transferred to five 100 µl drops of distilled water and letting it for 2 min on each drop. The sample was negative-stained with 1% uranyl acetate. The excess uranyl acetate was removed by contacting the grid edge with filter paper and the grid was air-dried. The grids were imaged with a JEOL 100CXII Transmission Electron Microscope operating at 100 kV. Images were captured on a Gatan Orius Digital Camera.
Nanoparticle tracking analysis (NTA). NTA was carried out using the ZetaView (Particle Metrix) following the manufacturer's instructions. EVs derived from saliva were further diluted 1000-to 5000-fold with PBS for the measurement of particle size and concentration.
Saliva proteomics (mass spectrometry). Saliva samples for proteomics processing were collected in a conical tube (same as EVs-see above) and stored in −80 • C on sample receipt. For proteomics processing, the Tandem Mass Tag (TMT) 6-plex kits were used (Thermo). Per sample, 300 µl of saliva were and dissolved in 300 µl lysis buffer (1:1 ratio saliva to lysis buffer to achieve > 2 mg/ml protein concentration). Protein concentration were evaluated using a Qbit (Life sciences). Per timepoint, 100 µg were used, adjusted to a final volume of 100 µl with 100 mM TEAB. The manufacturer's TMT labeling protocol was then followed to prepare the protein extract (part A, steps 7 onwards, and parts B for protein digestion and C for peptide labeling). Samples were ran through OffGel Fractionation (using an Agilent Offgel 3100 fractionator), and mass spectrometry was carried out with a ThermoFisher Q-Exactive mass spectrometer (www.therm o.com) using a FlexSpray nano-spray ion source. Survey scans were taken in the Orbi trap (70,000 resolution, determined at m/z 200) and the top twelve ions in each survey scan are then subjected to automatic higher energy collision induced dissociation (HCD) with fragment spectra acquired at 35,000 resolution. Additional details are provided in the online experimental protocols on Zenodo (see above).
Data mapping. Total saliva RNA-seq. Fastq files from paired-end sequencing (150 bp paired-end reads) were mapped using Kallisto 32,33 (with bootstrap sample parameter, -b, set to 100. For annotation, GENCODE 36 v28 transcripts and genome built GRCh38.p12 were used. The mapping results across timepoints were compiled using sleuth 34 (with DESeq 35 adjustment of Transcripts per Million). We note that the annotation used gene name concatenated with 'kind' information ('ext_gene':'kind').
The transcriptomics results were imported as OmicsObject constructs in MathIOmica 37 . Zero intensities were tagged as missing values, and intensities with aTPM < 1 were set to unity. Time series were constructed only for transcripts for which a signal was detected for at least 3/4 of the time points, and also constant time series were removed.
Total saliva proteomics. Proteomics .raw mass spectrometry files were analyzed using Proteome Discoverer (Thermo), using UniProt 38 human proteome database for reference. Mass tolerance was set to 10 ppm for precursor ions, and to 0.02 Dalton for fragment ions. Modifications included cysteine carbamidomethylation (fixed) and N-terminal and lysine TMT 6-plex and methionine oxidation (variable). Furthermore, we allowed for < 2 trypsin digestion missed cleavages. Proteins were identified using unique peptides of length ≥ 6 amino acids. We set FDR < 1% (strict) and < 5% (non-strict). For identification, we calculated results for both cases with 1 or 2 unique peptides per protein. We carried out peptide quantitation using unique peptides (reporter ion mass tolerance < 10 ppm ). For protein quantitation, we used medians of peptide ratios.
Multi-consensus reports from each set of technical replicates were constructed and used downstream in MathIOmica to construct an annotated OmicsObject. For each timepoint (sample) a Box-Cox power transformation was first used to transform the data to normal distributions 39 . Time series were constructed only for proteins with at least 2 unique peptides, and for which a signal was detected for at least 3/4 of the time points. Constant time series were also removed. www.nature.com/scientificreports/ EVs small RNA-seq. Small RNA-seq data from EVs were processed using the Genboree Workbench 41,42,83 exceRpt pipeline 40 to assess content by: (1) First removing reads that map to UniVec contaminants, 45S, 5S and mitochondrial rRNAs; (2) mapping reads sequentially to human miRNAs (mirBase), tRNAs (gtRNAdb), piRNAs (piRNABank), GENCODE and circRNAs (cirBase); (3) mapping unmapped reads from (2) to exogenous miRNAs and rRNAs; (4) finally mapping unmapped reads from (3) to all genomes in Ensembl and NCBI. Parameter settings and Genboree output files are available on Zenodo (see data availability). For each biotype MathIOmica OmicsObject constructs were created. Zero intensities were tagged as missing values, and intensities with aTPM < 1 were set to unity. Time series were constructed only for transcripts for which a signal was detected for at least 3/4 of the time points, and also constant time series were removed.
Temporal analysis and integration. For all mapped data time series were constructed with reference to the first timepoint for TFH1, TFH2 and TF , and with reference to the vaccination day for TFDaily. TF time series were constructed as paired hourly timepoint intensity differences between TFH2 and TFH1. All series were normalized as vectors to unit length. Time series classification analyses were carried out using MathI-Omica as detailed below and in the online Mathematica notebooks 37,84 .
Temporal classification details. The time series classification used MathIOmica's TimeSeriesClassification function with the Method -> "Autocorrelation" setting 37 . Briefly, for a given omic signal j with X j intensities over N times we construct a time series X j = X j (t 1 ), X j (t 2 ) . . . , X j (t N ) . The signal's periodogram is obtained using a Lomb-Scargle transformation [43][44][45] to account for uneven sampling, P LS . An inverse Fourier transform on P LS results in the autocorrelations ρ j for signal j as a list for lags 0 to n = ⌊N/2⌋ , ρ j = ρ j0 , ρ j1 , . . . , ρ jk , . . . , ρ jn . The autocorrelations' significance (p-value ≤ 0.01) is assessed by using a list of cutoffs, ρ c = ρ c1 , ρ c2, , . . . , ρ ck , determined from a null distribution of autocorrelations for each lag. These null distributions are generated from the calculated autocorrelations of simulated random signals that are created by bootstrapping (re-sampling of the original data with replacement). A signal is categorized in a class corresponding to the lowest lag deemed statistically significant, i.e. in class Lag l, where l = Min i : ρ ji ≥ ρ ci , and i ∈ 1, . . . , k . The result is a unique classification for each signal into Lag classes, which also ensures that any identified autocorrelation at a particular lag cannot possibly arise due to dependence on autocorrelations at lower lags.
The signals that do not show significant autocorrelation at any lag are checked for sudden signal spikes at any time point, and if so classified as spike maxima or minima. For each signal not showing autocorrelation, X j the signal maximum, max j = maxX j , and minimum, min j = minX j , are calculated across all time points. These values are compared against cutoffs {max cn , min cn } generated from bootstrap simulated distributions from the data. If for a signal X j of length n, max j > max cn it is classified in the SpikeMax class, or otherwise if min j < min cn it is classified in the SpikeMin class. Signals for which the signal intensity does not meet cutoff conditions are not reported.
Enrichment analysis. Gene-based over-representation analyses were run using MathIOmica 37 in Mathematica for GO and Reactome database entries. For miRNA enrichment analysis was run using the miRNA Enrichment Analysis and Annotation tool (miEAA) over-representation tool 85 .
Taxa groups were checked for over-representation using MicrobiomeAnalyst's 86 's web interface. Each subgroup and also each aggregate class were tested on multiple levels. The online MicrobiomeAnalyst database used included the following information, to show analyses at three different levels of mixed-level taxons, species-level and strain-level: • Mixed-Level Taxon sets included the following taxon sets: 1545 associated with host genetic variations, 239 associated with host-intrinsic factors such as diseases, 118 associated with host-extrinsic factors including diet and lifestyle, 446 associated with environmental factors such as drugs, chemical exposures and 53 associated with microbiome-intrinsic factors such as motility, shape, or spore forming. • Species-level taxon sets included: 61 associated with host-intrinsic factors including diseases, 92 associated with host-extrinsic factors including diet and lifestyle, 7 associated with environmental factors including drugs and chemical exposures. • Strain-level taxon sets included: 42 associated with host-intrinsic factors including diseases, 50 associated with microbiome-intrinsic factors such as microbe mobility and shape, and 399 associated with environmental factors including drugs and chemical exposures.
Network construction. Weighted expression networks were constructed in which each node represents one molecular species and each edge weight is defined as w ij = 1 (d ij +0.0001) , where d ij is the Euclidean distance between each pair of nodes {i, j} , and the offset 0.0001 was added to account for cases where d ij = 0 . Networks were constructed for both the classified TFD data and the TF data. To account for missing data in the computation of the Euclidean distance, mean imputation was used. Edge selection for the network construction was determined by filtering on one-tailed quantiles q(N) based on the w ij distribution in a given network k with N k nodes: www.nature.com/scientificreports/ Finally, in the network plots nodes were colored based on the MathIOmica classification group to which they belong.
Ethical approval. All experimental protocols were approved by the Institutional Review Board under protocol number LEGACY15-071 (15-071) at Michigan State University. All methods were carried out in accordance with the relevant guidelines and regulations. Informed consent was obtained from the participant as per the above protocol.

Data availability
Sequencing data reported are available on Gene Expression Omnibus under accessions GSE108664 (Saliva mRNA-sequencing) and GSE108666 (EV small RNA). Proteomics data are available on MassIVE as part of accession MSV000081869. All scripts and data analysis code utilized in the integrative analysis are available on Zenodo (DOI:10.5281/zenodo.3987587) as online data files (ODFs).
Received: 2 July 2020; Accepted: 23 December 2020 (1) q k (N k ) = 0. www.nature.com/scientificreports/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.