Volunteers’ microbiome had individual-specific composition

Twenty-six volunteers (13 women and 13 men) of age between 21 and 28 years participated in this study. On average 21 non-stimulated saliva samples (2 ml) were collected from each participant during the period of 30 days. DNA extracted from the saliva samples was amplified for the 16S ribosomal gene and sequenced. In total, 27,072,746 sequences from 551 samples collected from 26 volunteers have passed the quality filters and have been clustered into 18,146 operational taxonomic units (OTU). The final OTU names used in this study consisted from the assigned genus of the reference sequence and the number of the cluster.

The microbiome of all volunteers was dominated by the Streptococcus-OTU0 (highest similarity with Streptococcus parasanguinis) forming on average 37.5 ± 10.1% of the whole microbiome (Supplementary Figure S1) followed by taxon Rothia-OTU1 (20.9 ± 11.3%). The ratios of these two most dominant species differed significantly (p < 0.001, t-test) in 15% of pairwise comparisons among volunteers. For example, the proportion of Rothia-OTU1 could be as high as 35.3 ± 10.5% in volunteer F06 or as low as 5.4 ± 5.2% in volunteer M12. The overall microbiome composition of the volunteers differed significantly due to the less prevalent (average proportion 0.1–5.4%) individual-specific taxa (p < 0.001, “envfit” test, Supplementary Figure S2). Some volunteers (e.g., F18, F20, and M20) were characterized by high proportions of Granulicatella-OTU2 (8.7 ± 4.1%) and Atopobium-OTU3 (6.9 ± 3.7%), while other volunteers (e.g., M12, M21) had high proportions of Gemella-OTU6 (4.1 ± 3.1%). The volunteer F12 had the highest proportion of Saccharibacteria-OTU4 (8.9 ± 4.5%). Some volunteers did not possess any extreme proportions of the most detected prevalent OTUs (central part of the CCA plot in the Supplementary Figure S2). The indexes describing microbiome diversity, Shannon index (2.5 ± 0.4) and the evenness index (0.5 ± 0.1), also differed significantly (p < 0.001, t-test) in 11% of all pairwise volunteers combinations (Supplementary Figure S3).

Before commencing with the following analyses which may require larger computational memory, we used Procrust test to check whether reduction of the OTUs list would provide the same results as including all 18,146 detected OTUs (many of them may actually include OTUs corresponding to sequencing errors). The test showed that OTUs with proportions of above 0.1% (50 OTUs) result in the same ordination of samples as all detected OTUs (Supplementary Figure S4).

Volunteers differed by degree of microbiome temporal variability

Changes of the OTUs relative abundances over time was assessed by Taylor’s equation σ = V × μβ, where V and β are the parameters of the model, and σ and μ are the dispersion (in standard deviation) and the mean of the measurements. The variability parameter V was the direct estimator of the amplitude of fluctuations and therefore of the general stability, while the β parameter, the scaling index, described the statistical behavior of the ecosystem, always between 0.5 (behaving as a Poisson process) and 1 (behaving as an exponential distribution).29

All the values of the Taylor’s parameter β were between 0.5 and 1 which means that the most prevalent OTUs in all samples showed less relative variability over time than the less abundant OTUs (Fig. 1). The Taylor’s variability parameter V (average 0.27 ± 0.06) significantly varied in 26% of pairwise comparisons among volunteers (p < 0.001, t-test). The volunteer F11 had the least stable microbiome of our cohort (V = 0.47 ± 0.09), while the microbiome of F18 was the most stable (V = 0.20 ± 0.03).

Fig. 1 Temporal variability of the microbiome. Taylor’s Parameter space for the 26 volunteers of the study. V represents the y-intercept of the linear fit, and β to the slope of the line. Each individual has been placed in this plot according to its V and β value, where the error bars correspond to the standard error of the mean Full size image

In addition, we calculated the difference variability (DV) and rank variability (RV) to detect time-points with highest variability on intra-individual level (Supplementary Figure S5). DV expresses an absolute difference between every OTU’s rank (proportion) at a specific time point compared to the previous time point. The DV was generally higher in the least stable volunteer F11 when compared to the most stable volunteer F18 (Fig. 2). This is expected because there are more rank differences in subjects with higher Taylor’s parameter V (F11) than in subjects with lower fluctuations (F18). In general, DV and RV are good estimators of community shifts on intra-individual level. For example, the volunteer F04 had a short community composition shift during the days 20–24, but later on the day 25 came back to the previous bacterial composition (Supplementary Figure S5).

Fig. 2 Rank matrix for the 50 most abundant OTUs for F11 and F18. Rank matrix corresponding to the most and less time-variable volunteers, F11 and F18 respectively. Both plots represent the 50 most abundant OTUs, and heat-map colors corresponds to the abundance of each OTU at each time, ranging from light-yellow for rank 1, to black, representing very low ranks. Alongside, the Rank Stability Index (RSI) is colored by the percentage of rank stability, following the same color-code as before. Below, both Rank and Difference Variability (DV) is plotted in red and blue colors for each time point Full size image

The most prevalent OTUs Rothia-OTU1 and Streptococcus-OTU0 had the rank stability index (RSI) over 96% in both volunteers F11 and F18 (Fig. 2) which indicated that the most prevalent OTUs were very stable in the volunteer with the most stable microbiome as well as in the volunteer with the least stable microbiome. The microbiome temporal variability of the volunteers F11 and F18 differed due to the changing proportions of the less prevalent OTUs. The average RSI of all 50 OTUs in volunteers F11 (the least stable) and F18 (the most stable) were 52.3 ± 20.1% and 67.7 ± 19.2%, respectively. The volunteer F11 had 26 OTUs with RSI below 50%, while the volunteer F18 had only seven such highly unstable OTUs. When the OTUs were sorted according to their prevalence, the last OTU with an RSI above 70% in the least stable volunteer F11 was placed in the 16th position, while in the most stable volunteer F18 it was placed in the 49th position (Fig. 2). The OTUs found in lower proportions were in general less stable than the highly abundant OTUs (Supplementary Figure S5). However, the most prevalent OTUs were not always the most stable. For example, Saccharibacteria-OTU15 and Rothia-OTU15862 were in high proportion in F11 and F18, but it possessed a low RSI in the both volunteers (Fig. 2).

Large inter-individual and intra-individual difference of the oxidative stress markers levels

Oxidative stress can be measured by estimating oxidative damage to lipids (lipid peroxidation) and proteins (protein oxidation), or by quantifying the capacity to resist oxidative damage (antioxidant capacity).30 The lipid peroxidation was quantified by measuring thiobarbituric acid reacting substances (TBARS, average 0.10 ± 0.18 μmol/l, Supplementary Figure S6). Advanced glycation end products (AGEs, 0.27 ± 0.19 g/l) and advanced oxidation protein products (AOPP, 37.6 ± 21.8 μmol/l) were herein quantified as carbonyl and oxidative stress markers respectively, to express oxidative protein damage. The capacity to resist oxidative damage was measured by total antioxidant capacity (TAC, 578.6 ± 149.4 μmol/l) and ferric reducing ability of saliva (FRAS, 396.7 ± 183.6 μmol/l). The values of the five markers possessed intra-individual temporal variations, and they also differed significantly (p < 0.001, t-test) among volunteers—in 44–73% of volunteers pairwise combinations (Supplementary Figure S6).

In addition, the five oxidative stress markers were tested for their pairwise correlations on intra-individual level. FRAS correlated positively with TAC in 14 out of the 26 volunteers, but the majority of the salivary markers pairwise combinations resulted in non-uniform correlation patterns on intra-individual level among the 26 volunteers (Fig. 3). Correlations lagged by 1–3 days accounted for 36.7% of all significant correlations detected, while the remaining 63.3% were temporally direct correlations.

Fig. 3 Correlations between the oxidative stress markers. The summary of the Pearson’s intra-individual correlations. The color of square indicates either a positive or a negative correlation (blue or red) of a pair of oxidative stress markers in a volunteer. The number in square indicates the days by which the correlation was lagged. Empty white squares indicate no significant correlation. The majority of volunteers possessed a positive correlation between FRAS and TAC. Several marker pairs had positive correlations in some volunteers while negative correlations in others Full size image

Correlations between bacterial taxa and oxidative stress markers are individual-specific

As much as 94% of the 250 marker-OTU pairs (5 markers × 50 most prevalent OTUs) had a significant correlation on intra-individual level in some of the 26 volunteers (Fig. 4, details in the Supplementary Figure S7). These correlations were either temporally direct or lagged by 1–3 days. The remaining 6% of the marker-OTU pairs did not show any significant correlation in any of the volunteers. The possible delay in the correlations between our variables was measured for all possible combinations, and the significant results were only found on a single lagged day as Fig. 3 shows.

Fig. 4 Marker-OTU correlations on intra-individual level. The summary of the Pearson’s correlation plot in the Supplementary Figure S7. The graduated color in tones from blue to white to red indicate whether the correlations found in the volunteers on intra-individual levels were mostly positive or negative. The number on the left in each rectangle indicates the number of volunteers with a negative correlation, while the number on the right in each rectangle indicates the number of volunteers with a positive correlation. The oxidative stress markers and the OTUs have been ordered according to their correlation tendency by hierarchical clustering using Euclidean distances. The OTUs with the strongest tendency of either positive or negative correlations with some of the oxidative stress markers were e.g. Streptococcus-OTU16 (Streptococcus mutans), Actinomyces-OTU27 and Cardiobacterium-OTU41 Full size image

The results of the Pearson’s correlations showed that a particular marker-OTU pair can have positive correlations in some volunteers, while negative correlations in others (Fig. 4, details in the Supplementary Figure S7). This result is consistent with Spearman correlation coefficients, a non-parametric correlation (Supplementary Figure S7). The volunteers with contradictory correlation results did not possess any extreme values of a given marker-OTU pair.

Importantly, there were some marker-OTU pairs that produced mostly non-significant correlations, but when they were found to be significant in some of the volunteers, they were exclusively either positive or negative in these volunteers (Fig. 4, details in the Supplementary Figure S7). The most prominent exclusive correlations with oxidative stress markers had Streptococcus-OTU16 (S. mutans), a well established causative agent of dental caries.4 S. mutans correlated positively with AOPP (in two volunteers), with FRAS (in five volunteers) and with TBARS (in one volunteer) and correlated negatively with AGE (in three volunteers). The proportion of S. mutans in the salivary microbiome was as low as 0.19 ± 0.47% in the studied cohort indicating that its significant correlations with the oxidative stress markers were independent of its proportion in the microbiome.

Among other marker-OTU pairs with exclusively negative correlations found in a higher number of volunteers (3–6) were e.g., AOPP—Cardiobacterium-OTU41, AOPP—Actinomyces-OTU27, FRAS—Streptococcus-OTU56148 and FRAS—Cardiobacterium-OTU41. The marker-OTU pairs resulting in exclusive positive correlations were e.g. AOPP—Rothia-OTU15862, FRAS—Rothia-OTU18624, AGE—Streptococcus-OTU11787 and TAC—Actinomyces-OTU27 (Fig. 5).

Fig. 5 Marker-OTU correlations obtained by the simulation of collection of only one sample per volunteer. The bar-plots illustrate which portion (in %) of the 1000 combinations in the one-sample-one-volunteer approach simulation resulted in a positive correlation (blue), in a negative correlation (red) or no significant correlation (p > 0.05, grey). Both positive and negative correlations occurred in the majority of the marker-OTU pairs in low proportion. However, there were also some marker-OTU pairs with an explicitly positive or an explicitly negative correlation (no contradictory correlations were detected). The bacterial OTUs in this figure are arranged according to the Fig. 4, to make them easily comparable Full size image

Correlations and interactions between bacterial OTUs are also individual-specific

Similar to the marker-OTU pairs, the significant intra-individual correlations found among the OTU–OTU pairs did not show any notable generalized pattern (Supplementary Figure S7). There were some particular OTU–OTU pairs resulting in consistent correlations among most of the volunteers, for example a negative correlation was found between Streptococcus-OTU0 and Rothia-OTU1 in 17 out of 26 volunteers, while in the remaining 9 volunteers no significant correlation was detected.

Also, we analysed the interactions of the 15 most abundant OTUs from each volunteer using the generalized Lotka-Volterra model, a system of equations that has been often used for the inference of bacterial interactions in complex ecosystems such as the human microbiome.31 These models allow us, when a temporal series is available, to make better predictions of bacterial interactions rather than the use of classic correlation coefficients that could be filled by false positives and false negatives due to compositional effects.31 The results from the Lotka-Volterra model are more accurate than using simple correlations in terms of microbial interactions,32,33 but it did not yield any generalized interaction pattern either (Supplementary Figure S8).

Inconsistent results of the single sample collection simulation

In order to demonstrate the utility of the longitudinal sample collections in correlation studies, we simulated 1000 different single-samples combinations. Single samples from the 26 sample-sets were selected and combined to form a set of 26 samples in 1000 repetitions. Each of the 1000 simulated one-sample-one-volunteer sets was tested for correlations of oxidative stress markers with OTUs. The majority (82.4%) of the 250 marker—OTU pairs (5 markers x 50 OTUs) produced contradictory results. For example, a negative correlation was found in 22.2% of these single sample combinations for the TBARS—Rothia-OTU15862 pair, but still 0.1% of combinations for this marker—OTU pair resulted in a positive correlation, while 77.7% of the combinations did not yield any significant correlation (Fig. 5).

The results of the correlation analyses performed on this simulated one-sample-one-volunteer sets (Fig. 5) were not consistent with the results obtained on the intra-individual level (Fig. 4). For example, AGE—Actinomyces-OTU27 pair exhibited exclusively negative correlations in 26.6% of the one-sample-one-volunteer combinations, however, it showed positive correlations in two volunteers (F04, M01) on intra-individual level. Furthermore, the exclusive positive correlation of S. mutans (Streptococcus-OTU16) with FRAS observed in five volunteers on intra-individual level, was not confirmed by this one-sample-one-volunteer approach (only 2.2% of positive correlations and 0.3% of negative correlations in the 1000 combinations).

In contrast, in some cases the correlation results obtained by the two approaches were quite consistent; e.g., negative correlations between AGE—S. mutans, AOPP—Actinomyces-OTU27, AOPP—Cardiobacterium-OTU41, and FRAS—Cardiobacterium-OTU41.