Defining the resilience of the human salivary microbiota by a 520-day longitudinal study in a confined environment: the Mars500 mission

The human microbiota plays several roles in health and disease but is often difficult to determine which part is in intimate relationships with the host vs. the occasional presence. During the Mars500 mission, six crewmembers lived completely isolated from the outer world for 520 days following standardized diet regimes. The mission constitutes the first spaceflight simulation to Mars and was a unique experiment to determine, in a longitudinal study design, the composition and importance of the resident vs. a more variable microbiota—the fraction of the human microbiota that changes in time and according to environmental conditions—in humans. Here, we report the characterization of the salivary microbiota from 88 samples taken during and after Mars500 mission for a total of 720 days. Amplicon sequencing of the V3–V4 regions of 16S rRNA gene was performed, and results were analyzed monitoring the diversity of the microbiota while evaluating the effect of the three main variables present in the experimental system: time, diet, and individuality of each subject. Results showed statistically significant effects for either time, diet, and individuality of each subject. The main contribution came from the individuality of each subject, emphasizing salivary microbiota-personalized features, and an individual-based resilience of the microbiota. The uniqueness of Mars500 mission, allowed to dampen the effect of environmental variables on salivary microbiota, highlighting its pronounced personalization even after sharing the same physical space for more than a year. DRuCc4_oUAP4YMwxG5Lf3T Video abstract Video abstract


Introduction
The host-associated microbiota is stirring the attention from many fields of life science, including basic biology, evolutionary studies, biomedicine, and biotechnology. It is now well known that it plays several roles in modulating the host health and that changes in the composition of the microbiota in specific human body districts or organs (e.g., skin, gut, vagina, lung) may influence the correct functionality of other organs [1]. The concept of holobiont reflects the intimate relationships between the host and the microbiota [2,3], but it is often difficult to determine which part of the host-associated microbiota is in intimate relationships with the host vs. an occasional presence. Cross-sectional studies have been used to decipher the more stable, core, microbiota, present in all individuals analyzed, in comparison with the fraction which is more variable, i.e., present in few individuals only (see for instance [4]) and longitudinal analyses helped to understand the temporal stability of the microbiota [5].
The human microbiota is not a single entity, but it may have different characteristics and roles. The gut microbiota is expected to be more stable over time than other cavities that are more exposed to the environment-the oral cavity represents one of the first entry point of our body and is thus massively influenced by environmental conditions. The salivary microbiota is known to be affected by both biotic and abiotic factors [6,7], including the age, saliva chemical composition, tongue, and teeth [8]. Consequently, it is still under debate how much of the oral microbiota is stable over time and if this stability can be considered as a tight association with the host [7,9,10]. Given its sensibility to external perturbations, the salivary microbiota could be a good model to inspect the temporal dynamics and subject-by-subject variations impacting the human microbiota, but this sensibility could be a double-edged sword. Even if the disclosure of salivary microbiota temporal stability, and/or subject individuality, could indeed impact on scientific fields spanning from personalized medicine to forensic microbiology, controlling environmental exposures of salivary microbiota is difficult especially during our everyday life. Standardize these perturbations implies isolation procedures that are difficult to put in place.
Mars500 was the first long-term international study into interplanetary space flights. Managed by the European Space Agency and the Russian Space Agency, it was conducted in 2010-2011 when six male volunteers were kept for 520 days in a common confined environment established by the Institute of Biomedical Problems (IBMP) in Moscow, simulating a space flight to Mars. Data from Mars500 mission were studied from various point of view, including behavior [11], effect of cultural background [12], cognitive performances [13], circadian rhythms [14], hormone levels [15], and surface and gut microbiota [16,17]. Mars500 hence constitutes a unique experiment to determine, in a longitudinal study design, the composition and importance of the resident microbiota vs. a more variable microbiota (changing with time and environmental conditions) in humans. The aim of the work was to inspect the temporal dynamics of salivary microbiota, assessing the effect of diet regimes and individuality, using Mars500 as a unique long-term experiment where subjects were all confined in the same shared environment.

Salivary microbiota composition during the study
To inspect how the salivary microbiota reacts in a confined environment, we characterized samples collected during the entire duration of the Mars500 mission (720 days in total) by 16S rRNA gene amplicon sequencing of the variable region V3-V4 (Table 1). Table S1 summarizes all phases of the mission whereas Fig. 1a reports the sampling scheme used in this work. Amplified sequences formed 1890 amplicon sequence variants (ASVs) with a median number of 172.00 ASVs per sample (ranging from 81 to 317). A total of 4,337, 540 sequences specifically aligned to an ASV resulting in a sequencing depth ranging from 20,084 to 116,809 and a median value of 47,044 (for additional information about sequence analysis pipeline and the number of  Figure S3). Rarefaction curves reached a plateau above 15k reads suggesting an adequate sequencing depth for all samples (Supplementary Figure S4). Good's coverage estimator ranged between 99.99% and 100.00% across all samples indicating that roughly 0.01% of the reads in a sample came from ASVs that appear only once in that sample (Supplementary Table S5). Roughly, 99% of sequences aligned to variants that came from known bacterial taxa (Table S2). Supplementary Table S6 shows the overall taxonomic composition of samples whereas Figure S5 reports the phylogenetic tree reconstructed from ASVs. At phylum level Firmicutes, Bacteroidetes, Actinobacteria, Proteobacteria, and Fusobacteria accounted for more than 97% of the total number of reads assigned to taxonomically annotated ASVs (Fig. 1b and Supplementary Table S6). The total bacterial diversity (namely the alpha diversity) remained constant during the mission with no significant differences detected between the isolation period and the follow-up, across different diets, and across subjects (Table S7 and Fig. 1c-e). Also, time did not impact bacterial diversity as showed in Fig. 1f (random mixed model fitted using crewmembers as random intercept: slope lower than 0.002; Supplementary Table S8).

Effect of food and time on salivary microbiota
We inspected differences across samples (namely beta diversity) using non-metric multidimensional scaling (nMDS) on quantitative and qualitative indexes. Samples showed a similar distribution with all index tested (Fig. 2a): Sorensen index and unweighted UniFrac distance (qualitative analysis), and Bray-Curtis and weighted UniFrac distance (quantitative analysis). As opposed to alpha diversity, subjects, diets, and time significantly contributed to shape the salivary microbiota with different percentage of variance explained depending on the index but never exceeding 10% of the total variance (Fig. 2b). For all diversity indexes (except for the Sorensen index which reported a significant effect of subjects), the dispersion of tested factors was homogeneous, meaning that only the composition of samples varied among groups as highlighted by the permutational analysis of variance reported above. Diet impacted on bacterial genera usually present in the salivary microbiota of healthy subjects-such as Actinomyces, Veilonella, and Fusobacterium [18]-but also on Peptostreptococcus, Haemophilus, Megasphaera, and Prevotella, which have been correlated to different disorders of the oral cavity (such as periodontitis, dental caries, and oral lichen planus) [19][20][21]. Bacterial species classified as Alloprevotella, Fusobacterium, Dialister, Veilonella, and Megasphaera followed the same pattern: decreasing abundance passing from the first to the third diet and then back to starting values during the follow-up (even if the shift was not significant). Other species like Haemophilus and Prevotella reported a significantly higher abundance during a single diet-namely Haemophilus was more abundant during the follow up and the abundance of Prevotella was higher during the first diet. At phylum level, diets impacted more on Bacteroidetes and Firmicutes. Within them, 4 out of 28 (14%) and 3 out of 7 (43%) genera reported (at least) a significant difference during diet changes ( Figure S6).
To explore the effect of time on bacterial diversity, we used change-point analysis on both within-subject ( Figure S7 panel a) and between-subjects diversity ( Figure S7 panel b). Within-subject diversity measures changes in the salivary microbiota of each crewmember through time, whereas between-subjects diversity compares the salivary microbiota of different crewmembers at each time point (Fig. 3a, b). Three segments significantly divided within-subject diversity with two change-points at 123 days and 480 days. Between-subjects diversity was not segmented since the overall model gave better results than the segmented one according to the genetic algorithm used during optimization (Fig. 3b). The overall between-subjects model had an effect size of −0.00004 which means that after 520 days of isolation the overall diversity decreased by 0.02179. The effect of time on within-subject diversity was indeed higher than the (See figure on previous page.) Fig. 1 Salivary microbiota diversity along Mars500 mission and follow-up. a) Timeline of Mars500 mission. Phases were reported in the top part of the panel together with their length in days. The second phase included the landing simulation (20 days) as well as the trip back to Earth (250 days). Diets supplied to crewmembers were reported using different colors, violet for the first food variant, cyan for the second food variant, and yellow for the third food variant. Samples were reported using one point for each crewmember (subject 5001 red, 5002 pale blue, 5003 green, 5004 dark blue, 5005 ochre, and 5006 gray). The days of isolation were reported in the bottom part of the panel. b) Distribution of the main bacterial classes. Panels were divided according to crewmembers (vertically) and diets (horizontally). ASVs with a relative cumulative frequency lower than 5% in all samples were collapsed into a single group called "Other." c) Differences in alpha diversity-reported using the inverse Simpson index just like panels d, e, and f-between samples collected during and after the isolation period. d) Differences across diets (FV, first variant; TV, third variant; NR, normal diet). e) Differences among subjects. f) Differences along the whole mission and during the follow-up. Points are the average diversity values among subjects whereas error bars represent the 95% confidence interval around the mean. The red line represents the population effect of the linear mixed model whose coefficients are reported in Table S8 one observed for between-subjects diversity. During the first 123 days the effect modeled was 0.00103 reflecting an average increase of 0.12636 for all crewmembers. After the first change-point, within-subject diversity started to decrease with a regression parameter of −0.00064 (average decrease during the second segment of −0.22835). After the second change point, which roughly matched the end of the isolation period ( Fig. 3c, d), the within-subject diversity started to increase again. At the end of the follow-up period, diversity increased again of 0.29474 exceeding the average value detected in the first day of isolation ( Table 2).

Resilience of salivary microbiota
The average abundance of ASVs correlates with their persistence, the number of subjects in which a given ASV was detected at each time point. Figure   The dispersion of groups was tested for homogeneity and results were reported on the top of each ordination (a pvalue higher than 0.05 means that dispersions are homogeneous). b) Permutational multivariate analysis of variance using distance matrices based on the same indexes reported in panel a. The R 2 values associated with each factor used in the analysis is reported in the horizontal axis whereas asterisks report the significance level of each factor (*, p-value < 0.05; **, pvalue < 0.01). Colors represent the different factors modeled in the analysis. For additional information about diets, sampling point and crewmembers see Supplementary information and Supplementary Figure S1 and Figure S8, panel a and b). The inconsistent microbiota showed low average persistence in respect with the stable microbiota, but it contained the largest amount of variants (1746 ASVs against 144 of stable microbiota). Unlike stable ASVs, subjects lost and acquired inconsistent ASVs both during and after the isolation period ( Figure S8 panel b and c). Stable ASVs were detected in roughly 30% of all subjects at each time point (26 samples on 88) with sporadic losses and acquisitions ( Fig. 4a and Figure S9).
We represented the acquisition and loss of bacterial species during the whole mission using networks. At each time point, we linked subjects to ASVs detected in their salivary microbiota forming a bipartite network structure which reflected the underlying bacterial community structure. The loss and acquisition of bacterial ASVs were shown in supplementary video S1 where green squares represent subjects, red circles represent inconsistent microbiota, and light blue circles represent stable microbiota. As shown in the video, the topology of the networks did not change in time, but at each time Crewmembers' salivary microbiota composition in time. a) Bray-Curtis (also known as quantitative Sorensen) index has been used to inspect distances between and within-subject during isolation and follow-up. Between-subjects diversity was computed by comparing the salivary microbiota of each subjects at each timepoint (gray arrows); within-subject diversity was computed by comparing the salivary microbiota of the same subject over time (black arrows). b) Change-point analysis revealed changes in salivary microbiota composition of each subject (CP detection). Genetic algorithm and linear modelling detected increasing/decreasing patterns along time (GA optimization). Finally, we fit a linear mixed-model for each segment detected using crewmembers as random intercept (Modelling). c and d) Results obtained following the pipeline reported in "b" for within-and between-subjects differences. Diets were reported in the bottom part of the plots using different colors (FV, fist food variant; SV, second food variant; TV, third food variant). Since crewmembers ate freely during the follow-up, no diet was reported in the plot Table 2 Temporal changes of salivary microbiota. Within-subject diversity was divided into three segments following change-point analysis whereas between-subjects diversity was modeled on the full-time period since no change-points were detected. Results of mixed effect models fitted for each segment were reported in the table. b, regression parameter (slope of the model); SE, standard error; t, t-value (also known as "standardized" regression parameter); df, degrees of freedom; p, p-value.  Table S9). The end of the isolation period significantly increased the average number of formed edges (namely acquired ASVs) of 28 but the trend was still negative (Fig. 4b). The number of formed edges was independent from the number of lost edges (Spearman's ρ = −0.04; p value = 0.750). The salivary microbiota structure did not change at the beginning, during, and after the isolation period reporting a similar network topology ( Fig. 4c and supplementary video S1). Amplicon sequence variants clustered into cluster 1 had a marginal position in all networks, linking only one (or a few) subjects-in many cases no link was reported for ASVs of cluster 1 since they are intermittently present in most subjects as reported in supplementary video S1. In contrast, ASVs clustered into cluster 2 had a more central position in all community networks reporting connections with at least four subjects in the majority of the cases. The topology of the networks was confirmed by the centrality analysis performed on ASVs of both clusters 1 and 2. The centrality of ASVs assigned to cluster 2 was higher than those assigned to cluster 1 at every time point, highlighting the central role of these variants in respect with the whole community structure (Fig. 4d, Wilcoxon rank sum test: all p-values < 0.01). Also, the abundance reported a similar effect with ASVs of cluster 2 showing higher values during the whole experiment (Fig. 4e, Wilcoxon rank sum test: all p-values < 0.01).

Drivers of diversity
To inspect drivers of beta diversity along and after the isolation period, we fitted a linear model for each ASV detected in the salivary microbiota of crewmembers. The Mars500 mission time-scale was divided into three stages according to the changepoints detected for within-subject diversity. Crewmembers showed a different number of ASVs reporting a trend of diversity similar to the one reported in Fig. 3c. The number of ASVs showing a significant effect of time and changepoints ranged from 3 to 28 depending on the subject ( Fig. 5a and Table S10). Stable ASVs-namely those grouped into Cluster 2 according to their prevalence in timethat showed a significant trend of diversity enriched the saliva of four out of six crewmembers, if compared with the overall occurrence (Fig. 5b). The fraction of stable ASVs in subjects 5004 and 5005 was roughly six times higher than the average fraction of stable ASVs, whereas  Fig. 3c. Lines show the predicted differences using linear modeling. We standardize differences to represent ASVs with different ranges of values in the same panel. b) Enrichment analysis of members of cluster 2 in respect with the population. The firs bar on the left represents the overall fraction of AVSs assigned to cluster 2 whereas the other bars report the fraction of ASVs assigned to the same cluster for each subject. Adjusted p-values were reported using a single asterisk (p < 0.05) or two (p < 0.01) subjects 5002 and 5006 reported a fraction three times higher than the average.

Discussion
The reported experimental results under Mars500 mission evaluated the temporal dynamics of the human salivary microbiota in a controlled and confined environment. All samples from the 6 (male) crewmembers included Firmicutes, Bacteroidetes, Actinobacteria, Proteobacteria, and Fusobacteria as main phyla and conserved their taxonomic composition along time and among individuals, as previously reported by other studies (see for instance [7,9,10,[22][23][24]). Despite this conservation, external factors significantly influenced the salivary microbiota composition, but their influence was restricted to a low percentage of the community (less than 10% of total variance explained). We found a great number of species that intermittently passed through the salivary microbiota but they neither affected its overall structure nor its taxonomic composition. Despite their number, these transient species struggle to thrive in human saliva reporting a low abundance during the whole experiment. The presence of these species was limited to one or a few crewmembers in each time point and their centrality was extremely low, suggesting a marginal role in maintaining the whole community structure-that was highly conserved throughout the experiment as witnessed by slight changes in the topology of the community networks at different time points. Inspecting the oral microbiome of completely isolated subjects in time helped to discern between undersampled ASVsbacterial species that fall under the detection threshold and are thus not detected by sequencing methods even if they are present in the oral cavity-from bacterial species that are truly absent from the oral cavity of crewmembers. Mixed effect models helped to reduce the error and detect transient ASVs with a higher accuracy but discerning from udersampled and really transient species remains a complex task that needs to be confirmed by other molecular methods. On the other hand, a hundred bacterial taxa dominated the salivary microbiota composition of all crewmembers, rarely changing host even when the isolation period ended. The isolation time affected bacterial diversity of single individuals, but it did not alter the microbiota of the whole crew. The bacterial diversity of crewmembers decreased-if compared within consecutive time points of the same individual-but the effect ended immediately after the isolation period when it started to increase again. External perturbations, impossible to control outside the isolation facility, modulated the salivary microbiota composition when crewmembers got out the isolation facility and likely started eating different types of food, getting in touch with other people, or simply visiting different places.
The findings obtained under Mars500 mission suggest that sharing the same confined environment-and possibly following the same diet regime with few variations-imbalanced the ecology of the salivary microbiota while reducing its complexity. Unfortunately, we do not know if this effect could-positively or negatively-affect the health of the hosts. Extending this concept, we could suggest that a depauperate microbiota is less reactive to sudden changes in external conditions weaken the host adaptive capacity. However, a loss of strains in the microbiota corresponds to a loss of putative members that could be relevant under changing environmental conditions [3,4,25]. The diet regimen used during the study-although it coincided with the isolation period, and thus, the two effects could not be mutually excluded-played a role in the loss of species and complexity reduction but the intra-subjects variability had a higher impact in shaping oral the oral community. The isolation reduced crewmember oral microbiome complexity on an individual basis, highlighting that a shared diet regimen may reduce the diversity of a single subject in time, but it could poorly affect the whole crew. Between-subjects differences remained unaltered suggesting that individuals follow somewhat independent dynamics of their salivary microbiota (i.e., personalized dynamic). The same evidence indicates that crewmembers, though sharing the same environment, did not exchange their salivary microbiota, leading to hypothesize that quite stable personal salivary microbiota features are present in humans and confirming the betweensubject effect. The presence of sampling points before the isolation period could have enhanced our understanding of these personalized dynamics by setting up a "baseline microbiota" for each crewmember. Unfortunately, as for the other experiments done within Mars500 mission (see for instance [16,17,[26][27][28]) due to the military nature of the Mars500 experiment, no sample of crewmembers before the mission was taken/ made available, limiting our discussions to the isolation period itself. Other studies based on the Mars500 mission had the same experimental design but came to slightly different conclusions probably due to different human districts sampled. The fecal microbiota of the same crewmembers showed an increasing trend of similarity among subjects, especially in relation to a sharing of rare taxa [17], indicating that salivary microbiota has a more pronounced personalization than fecal microbiota. Indeed, from our data, there is no evidence of an ecological succession, as those shown in the fecal microbiota of the same crewmembers. In support of a more pronounced personalized dynamics of salivary microbiota than of fecal microbiota in Mars500 mission is the finding that the ASVs contributing to the temporal dynamics in the subjects were highly variable (from 3 to 28) indicating that each crewmember salivary microbiota changes was driven by a peculiar (viz. personal) set of microbial taxa.
Even if factors such as diet and time influenced the salivary microbiota composition of crewmembers, most details of the Mars500 missions are unknown. The mission was a military experiment and several outcomes are still sealed. The composition of diets for example is unknown, and thus, any speculation on the effect of particular food intake would not be grounded. Despite these limitations, the use of 16S rRNA gene amplicon sequencing, allowed us to detect key features of human salivary microbiota under a condition that is almost impossible to replicate. The complete isolation of the participants of the mission made possible the first observation of salivary microbiota composition minimizing the effect of external perturbations.

Conclusions
The reported longitudinal analysis of human salivary microbiota confirmed the stability of the microbiota over time and suggested the presence of resilient personalized taxonomic features, which may deserve further attention. This study allowed to elucidate the contribution of a stable and confined environment, as that of Mars500 mission, in reducing the microbiota diversity while controlling for the effect of diet on salivary microbiota.

Methods
The Mars500 experiment Mars500 mission was conducted in 2010-2011 by Russia's Institute of Biomedical Problems (IBMP), with extensive participation by the European Space Agency (ESA) as part of the European Programme for Life and Physical Sciences (ELIPS) to prepare for future human missions to the Moon and Mars. The whole project consisted of three isolation studies: a 14-day pilot study to test facilities and procedures used during the simulation, a 105-day pilot study involving six crewmembers, and a 520-day study that simulated a complete space flight to Mars and back. Mars500 crew was composed of six male volunteers. All crewmembers were confined in the same living space from the 3rd of June 2010 till the 4 of November 2011 when they finally stepped out of the isolation facility to come back to their normal activities. During the mission, the crew was hermetically isolated from the rest of the IBMP facility. Crewmember received three type of diets, a so called "first variant" (FV) and "third variant" (TV) and, after the experiment, returned to a normal (non-supervised) diet regime (NR). Three crewmembers, which participated to the landing simulation on mars, received another type of diet (called "second variant"), but no sample were collected during the mission due to the absence of sampling facilities in the landing module (from the 251st till the 270th day of the whole mission). Detailed information on the experiment is reported in Supplemental Information file. Further details are also reported in the companion Mars500 microbiology paper [17].

Collecting salivary samples
Saliva samples were collected individually, based on the scientific protocol and pre-confinement training, with 5ml sterile vials (Nalgene V5257-250EA). Upon completion of the saliva sampling, all samples from one sampling event were put into the hatch. After that, they were removed by the responsible person of the IBMP and stored at −80°C. After being stored at −80°C in the laboratories of the IBMP for periods of at least 4 days up to 6 months, the samples were sent via World Courier. Shipping from Moscow to the University of Tuscia, Viterbo -Italy. The shipment was performed in three batches on dry ice to avoid repeated freeze-thaw cycles which lead to reduction of microbial viability. Upon arrival, samples were kept at −80°C until processing. Salivary samples were collected by crewmembers during and after the permanence in the isolation facility (Table 1): 42 samples were collected during the first simulated journey from Earth to Mars (seven time-points), 30 samples were collected during the simulated trip back home from Mars to Earth (five timepoints), and other 16 samples were collected when crewmembers came back to their normal activities and were followed for additional 200 days (three time-points). Unfortunately, two samples collected in the latter stage gave no good quality DNA and were thus discharged. Crewmembers did not collect salivary samples during the simulated landing on Mars where three of them-which simulated the landing on a separate module-used a different food variant. For this reason, the second food variant used was not reported in the work passing from the first food variant (FV) directly to the third food variant (TV).

Amplicon sequence variant reconstruction
The DADA2 pipeline (version 1.14.1) was used to reconstruct amplicon sequence variants (ASVs) from Illumina reads [32]. Both ASV reconstruction and statistical analyses were performed in the R environment version 3.4.3 (http://www.R-project.org). For a complete description of all the step performed see Supplementary material section. Briefly, primers used for V3-V4 amplifications were detected and removed using cutadapt version 1.15 [33]. Low-quality reads were discarded using the "filter-AndTrim" function with an expected error threshold of two for both forward and reverse read pairs. Denoising was performed using the dada function after error rate modelling ("learnErrors" function). Denoised reads were then merged discarding those with any mismatches and/ or an overlap length shorter than 20bp ("mergePairs" function). Chimeric sequences were removed using the "removeBimeraDenovo" function. Taxonomical classification was performed using the IDTAXA algorithm [34] available in the DECIPHER package version 2.14.0 against the latest version of the pre-formatted Silva small-subunit reference database (SSU version 132 available at http://www2.decipher.codes/Downloads. html) [35,36]. All sequences classified as chloroplasts, mitochondria, Archaea, and Eukarya were removed. A summary of retained reads in each step is reported in Table S1 and in Figure S1.

Diversity estimation
Bacterial diversity in each sample was computed using inverse Simpson index as implemented in the "diversity" function of vegan package. Differences according to crewmembers, permanence in the isolation facility, and food variants were inspected using one-way analysis of variance (ANOVA). The effect of time was modeled using linear mixed models with fixed slope and random intercept. Since alpha diversity was measured multiple times on the same statistical units, crewmembers were used as random intercept factor.
Diversity across samples was inspected using different approaches. Qualitative and quantitative indexes were used to infer pairwise distances between samples. Qualitative indexes are binary indexes which take into account presence/absence of species to compute distances between samples whereas quantitative indexes are mainly based on the abundance of species [37]. Sorensen index [38] and un-weighted UniFrac distance [39] were used as qualitative indexes whereas Bray-Curtis dissimilarity [40] and weighted UniFrac distance [39] were used as quantitative indexes. UniFrac distances were computed using the distance function of the phyloseq R package version 1.30.0 [41] whereas Sorensen and Bray-Curtis dissimilarity indexes were computed using the "vegdist" function of the R package vegan version 2.5-6 [42]. Differences between salivary microbiota composition of the same crewmember at consecutive time-points (within-subject diversity) were computed using the "TBI" function of the adespatial R package version 0.3-8 [43]. Packages vegan and adespatial use the same definition of Sorensen and Bray-Curtis distance.
Sorensen index is defined as follows: where A and B are the numbers of ASVs on compared samples, and J is the number of the ASVs shared by both samples.
Bray-Curtis index is defined as follows: where b ik is the Bray-Curtis distance between sample i and k with J number of species, n ij is the abundance of species j in sample i, and n kj is the abundance of species j in sample k [44]. The adespatial package computes the index by splitting the total diversity into 3 components: A (the unscaled similarity between two samples), B (the species loss between two samples), and C (the species gain between two samples). The computation is then performed using the formula (B + C)/(2A + B + C) giving the same results as the formalization given above. Qualitative indexes rely on the assumption that all taxa are equally contributing to bacterial diversity independently from their abundance. For this reason, even extremely rare taxa may be relevant in shaping sample distribution. To relax this assumption ASVs detected in less than 5% of samples (4 samples overall) with an abundance lower than 10 were filtered out before diversity calculation.
Distances across samples were reported using nonmetric multidimensional scaling (nMDS) as implemented in the "metaMDS" function of the vegan package, with 300 random starts and monotone regression [45,46]. To test the effect of food variants, crewmembers, and time in shaping the salivary microbiota, we used permutational multivariate analysis of variance on distance matrices obtained above ("adonis2" function of the vegan package with 1000 permutations). The proportion of sum of squares from the total (namely the R 2 value of permutational analysis) was used to report the percentage of variance explained by each factor included in the analysis. Before testing for differences in bacterial composition among groups is advisable to make sure that groups are homogeneously dispersed; otherwise, permutational tests (such as adonis) may report significant results entirely due to uneven dispersion. To distinguish between actual differences in composition or differences due to dispersion, we used the "betadisper" and "anova" functions (vegan package). P-values obtained were corrected using Benjamini & Hochberg correction (also known as false discovery rate) [47]. To avoid possible biases induced by uneven sequencing depths, read counts were scaled using DESeq2 before diversity calculation ("counts" function) [48]. Scaled counts were additionally transformed using the square root of the Wisconsin double standardized counts ("wisconsin" function of the vegan package).

Influence of time on bacterial diversity
Salivary microbiota may be affected by several factors. Sharing the same environment for a prolonged period of time may alter the composition of salivary microbiota at different levels. The bacterial composition may be altered within the same individual taken at consecutive time points but even between multiple individuals at each time point. To inspect both of these components, bacterial diversity within and between subjects, calculated as reported above, was modeled through time. Change-point analysis was used to identify specific time points which led to a decrease/ increase of diversity. The optimal positioning and number of change-points for each crewmember was identified using a non-parametric cost function as implemented in the "cpt.np" function of the changepoint.np R package, version 1.0.2 [49]. The pruned exact linear time algorithm (PELT) [50] was used to detect temporal changes in diversity within the same subject and between different subjects. The PELT algorithm searches for an optimal solution by minimizing the cost of different segmentation. We used the modified Bayes information criterion penalty term (MBIC) [51] as penalty function for cost minimization by the algorithm. Since PELT algorithm is exact, a solution is always found for each time series so, to avoid inflation of change-points due to the presence of data coming from six different subjects, a genetic algorithm was used to fine-tune the analysis. All change-points detected were used as starting point of the genetic algorithm and a fitness function was defined as: where, RMSE(lmm k ) stands form the root mean square error of the generalized linear model constructed on the k segment and n is the number of segments defined by change-point analysis. High error corresponds to a low fitness value whereas low error corresponds to a high fitness value. At each step of iteration, the genetic algorithm will keep segmentation that led to linear models with a low RMSE and discard those leading to high error models. We implemented this algorithm using the R package GA [52] version 3.2 with a population size of 200. Generalized linear model were fitted using crewmembers as random intercept and pvalue were computed with the Satterthwaite's degrees of freedom method as implemented in the package lmerTest [53,54] version 3.1-2.
We assessed the persistence of the salivary microbiota across subjects using dynamic time warping algorithm implemented in the dtwclust R package (version 5.5.6) [55]. At each time point, the number of subjects in which a given ASV was detected (namely the ASV's persistence) was reported together with its abundance. Persistence matrix was scaled and centered before clustering. Centering was performed by subtracting the mean of each ASV from their persistence whereas scaling was performed by dividing each value by its standard deviation. The relation between persistence and abundance was tested by fitting a linear model ("lm" function of R stat package). All details about clustering and modeling were reported in Supplementary methods. Persistence across subjects was also used for network construction: at each time point, we constructed a bipartite network by linking subjects with ASVs that were present in their salivary microbiota. Subjects were represented using squared nodes whereas ASVs were represented using round circles colored according to the groups defined above. Doing so, we generated twelve bipartite, acyclic, and undirected networks representing the salivary microbiota of all subjects at different time points. The network R package (version 1.16.0) was used for network reconstruction whereas the package ggnetwork (version 0.5.8) was used for plotting [56,57]. Node-level centrality score was computed using the "centr_degree" function of the R package igraph (1.2.6). The effect of time and the end of the isolation period on the number of formed/destroyed edges were tested using mixed effect models with random intercept. Subjects were taken as random intercept whereas the time and the end of the isolation period as fixed effects. P-values were computed as discussed above.
Differences along time, for each ASVs detected, were inspected by selecting drivers of within-subject beta diversity. Absolute differences between consecutive time points were fitted using linear models and the effect of time and changepoints was inspected. The slope of the models was used to assess the trend of selected ASVs in each stage. We selected ASVs showing a trend similar to what we found during changepoint detection to focus the analysis only on components of the salivary microbiota of each subject contributing to the total diversity. Finally, to inspect if "drivers of diversity" detected were more present in stable or inconsistent salivary microbiota, a hypergeometric test was performed. The test compared the occurrence of stable ASVs (Cluster 2) in the overall population against the microbiota of single