Gut microbiota ecology: Biodiversity estimated from hybrid neutral-niche model increases with health status and aging

Volkov treats this model in the framework of the Chemical Master Equation and proves that under these assumptions the probability distribution of the RSA is expected to be the zero-truncated Negative Binomial PRSA n ð Þ ¼ Nobs 1 ð1 b=dÞS=b ð1 b=dÞS=b GðS=bÞ ðb=dÞn n! G n þ S=b ð Þ ð8Þ where Nobs refers to the total number of observed species and Γ is the gamma function. Biodi- versity can be estimated using the Hubbell biodiversity index [30, 31], defined as y ¼ Nobs ½ð1 b=dÞ S=b 1� � GðS=bÞ ð9Þ We modelled the empirical RSA derived from 16S rRNA data considering three possible scenarios. First, we tested pure neutrality by fitting the data with the 1NB model (Eq (8)). Then, we relaxed the hypothesis of species equivalence considering a hybrid niche-neutral model (2NB model) that assumes the existence of two non-interacting neutral niches (the evo- lutionary dynamics of each niche is neutral). Finally, we further relaxed the neutral hypothesis contemplating a hybrid niche-neutral model with three niches (3NB model). The mathematical aspects of the 2NB model are detailed in Bazzani et al. [38] and show that the stationary RSA distribution is a mixture of two zero-truncated Negative Binomials PRSA n ð Þ ¼ a � y1 ðb1=d1Þn n! G n þ S1=b1 ð Þ þ 1 a ð Þ � y2 ðb2=d2Þn n! G n þ S2=b2 ð Þ ð10Þ where α is the mixture coefficient, bi, di and Si are the birth, death, and influx rates of the i-th niche (i = 1,2), and θi is the biodiversity number relative to niche i, and is equivalent to the one in Eq (9).

Analogously, we assume the stationary RSA distribution for the 3NB model to be a mixture of three zero-truncated Negative Binomials (Eq (11)) to which 3 biodiversity numbers are associated: θ1, θ2 and θ3. PRSA n ð Þ ¼ a � y1 ðb1=d1Þn n! G n þ S1=b1 ð Þ þ b � y2 ðb2=d2Þn n! G n þ S2=b2 ð Þ þ 1 a b ð Þ � y3 ðb3=d3Þn n! G n þ S3=b3 ð Þ ð11Þ Model fitting and model selection Empirical RSA distributions were computed by counting the number of OTUs with a certain number of individuals. The 1NB, 2NB and 3NB models were fitted to the data using a custom implementation of the Approximate Bayesian Computation (ABC) rejection algorithm. Assuming a degree of similarity across samples, the goal was to implement a hierarchical model to estimate the parameters of the model, allowing the estimate from the general popula- tion to inform about the plausible values of each individual sample.To facilitate the algorithm convergence, we implemented this method with a two-step approach. In the first step, we assigned uninformative prior distributions to the model parameters and we fitted the model for each sample using the ABC algorithm detailed below. Then, for each data set, we con- structed the posterior distributions of the model parameters pooling the accepted parameters of all samples belonging to that data set and fitting them with either a Beta distribution (mix- ture coefficients) or a Gamma distribution (all the other parameters). Notice that, in general, PLOS ONE Gut microbiota biodiversity and healthy aging

PLOS ONE | https://doi.org/10.1371/journal.pone.0237207 October 30, 2020 15 / 23

the number of accepted parameters was highly variable from sample to sample. For some sam- ples the number of parameters accepted in the first step was very small or null. For other sam- ples, however, thousands of parameters were identified as acceptable. For this reason, we chose to limit the number of accepted parameters obtained from each sample to 5. This allowed to avoid high unbalances in the number of parameters derived from each sample, as well as to improve the computational efficency of the posterior fitting. Finally, the posterior distributions obtained in the first step were used as data set-specific prior distributions and the ABC algo- rithm was run again to obtain the final model estimates (second step). In both fitting steps, for 107 times we randomly sampled a set of parameters from their prior distributions and we sim- ulated a number of data equal to the number of OTUs according to the selected model (1NB, 2NB or 3NB). Then, we evaluated the accordance between the simulated and the empirical RSA comparing the two Preston’s plots. This choice was made to reduce the numerical issues related to the sparsity of the data in the heavy tailed distributions [30]. Specifically, each set of parameters was accepted if both the conditions described in the following were satisfied. First, we set the maximum acceptable absolute difference between data and simulation counts in each bin to 30% of the data counts. Secondly, we constructed a variant of the chi-squared test that was appropriate for the comparison of the two observed samples (real and simulated

data). Our null hypothesis, here, is that the difference between the observed values for each bin of the Preston’s plot is distributed as a Skellam’s distribution (the difference between two sta- tistically independent random variables, each Poisson-distributed with respective expected val- ues μ1 and μ2). Given that this distribution converges to the Normal distribution, we can use a variant of the chi-squared test obtained by summing the standardized Skellam’s distributions of each bin of the Preston’s plot: w2 ¼ PN i¼1 ðmi1 mi2Þ2 ðmi1 þ mi2Þ ð12Þ where i indicates the i-th bin of the Preston’s plot. The number of degrees of freedom was set equal to the number of non-zero bins and the set of simulated parameters was accepted when the probability that the simulated and true bins of the Preston’s plot come from the same distribution was higher than the probability that they did not, i.e. when the chi-squared cumulative distribution at χ2 was lower than 0.5. This crite- ria was designed according the samples acceptance criteria in ABC. Finally, for each sample, the posterior distributions of the parameters were computed considering all accepted set of parameters. Model selection was achieved comparing the posterior probability of each model, estimated computing the ratio between the number of accepted parameters sets over the total number of simulated parameters (107), obtaining a median number of accepted parameters of ~2000 sam- ples (acceptance rate of ~2�10−4, relative error in the rate estimation of 2%). Specifically, for

each pair of models we computed the logarithm of the ratio of the two posterior probabilities of each sample. This method computes an approximation of the Logarithmic Bayes Factor, defined as the log ratio of the probabilities that each model is the actual true model, without incorporating prior beliefs about the plausibility of each mode [47]. This method is akin to a Bayesian equivalent of the likelihood ratio tests such as BIC (Bayes- ian Information Criterion) and AIC (Akaike Information Criterion) [48], but including an implicit penalization for the number of parameters that each model possesses and the shape of the prior for each parameter: wider priors (encoding for less certain parameters) cause a greater penalization than narrower ones (encoding for more information available about them). PLOS ONE Gut microbiota biodiversity and healthy aging PLOS ONE | https://doi.org/10.1371/journal.pone.0237207 October 30, 2020 16 / 23

All the methods were implemented in python 3.6.8 [49] using the following libraries: patsy [50], pandas [51], numpy [52], scipy [53], matplotlib [54], seaborn [55], pymc [56]. The code to perform the parameter estimation through ABC starting from the OTU abundances table is available on a github repository at https://github.com/UniboDIFABiophysics/ AlphaDiversityPublic. Classical biodiversity indices Results obtained with θrare and θabundant were compared to those obtained when biodiversity was computed using classical biodiversity indices that are based on the empirical relative abun- dance of OTUs. Specifically, we considered Shannon [25], Pielou [27], Hill1 [28], Simpson [26] and Hill2 [28] indices, whose definitions are reported in the following. Let us call pi the proportion of individuals belonging to the i-th OTUs and S the total num- ber of OTUs. Then, 1. Shannon index (entropy) is defined as H = −∑pi�ln(pi) 2. Pielou index (eveness) is defined as J = H/Hmax, were Hmax is the maximum possible value of Shannon index H, that is achieved when pi = 1/S for each i, i.e. Hmax ¼ P 1 S � ln 1 S � 3. Hill1 number (Hill number of order 1) is defined as Hill1 = exp(H) 4. Simspon index is defined as l ¼ P p2 i . Notice that differently from the other indices Simp- son index decreases with biodiversity. 5. Hill2 number (Hill number of order 2) is defined as Hill2 ¼ 1 l Spline regression: Modeling the trend of GM biodiversity with aging The general relationship between the GM biodiversity indices, θrare and θabundant, and aging

was investigated by fitting a natural cubic splines model [57] adjusted for sex (0 = male, 1 = female) and standardized total number of reads, using the library patsy [50] in python 3.6.8 [49]. See the following paragraph “Notes on splines regression” for details. Here, we considered the healthy control subjects from all the data sets except the Italian controls of Schnorr et al., for which age was not available, and the ELDERMET data sets at times greater than 0, for which young controls were not available. Samples with outlier total number of reads (z-score > 4) were removed prior the analysis. These included 6 elderly sub- jects from the ELDERMET study at T0 data set, 8 subjects from the Bian data set, 4 subjects from the Odamaki data set, and 1 subject from the Kong data set. The GM biodiversity estimates, θrare and θabundant, were analysed separately. After standard- izing the diversity indices within each data set, we fitted the natural cubic spline model consid- ering the data set to which the samples belong as confounding variable. The number of degrees of freedom was set to 3 based on the Akaike Information Criteria (AIC) and the Bayesian Information Criteria (BIC), as shown in S3 Fig. The spline regression model was also computed considering classical biodiversity indices, i.e. the Pielou, Shannon, Simpson, Hill1, and Hill2 indices. Notes on splines regression Splines regression models (also referred to as Generalized Additive Models [57]) allow to describe non-linear behaviors in the data using the framework of linear regression, where the

outcome variable is described as the linear combination of generic, non linear functions of the independent variables. While traditionally this was performed combining functions such as PLOS ONE Gut microbiota biodiversity and healthy aging PLOS ONE | https://doi.org/10.1371/journal.pone.0237207 October 30, 2020 17 / 23

squares and cubes of the data, the spline regression uses a different set of basis function, one of which is always the average value of the outcome variable, usually referred to as the intercept, albeit improperly. The other functions are chosen from a family of functions depending on the desidered properties. In the case of biological functions, a common approach is to fit with a polynomial (in this case of 3rd order) a subset of the data, while imposing 4 conditions: continu- ity of value, continuity of derivative of first order, continuity of the derivative of second order, and the overall function at the boundaries is linear. These are usually referred to as natural cubic splines. The number of knots (the junction points between different parts of the data set) were chosen using a cross-validation approach, using the AIC criterion as the chosen metric. Due to the constrains required for this functions to create a smooth fit, they can exhibit small fluctua- tions outside the original domain, such as those that can be seen in spline 1 and 3 in Fig 1. Classification and prediction of health status in old age Elderly samples from the ELDERMET data set were classified as healthy or unhealthy accord- ing to the FIM score [58], Barthel index [59] and MMSE score [60] and to the residence set- ting. In particular, for each subject, a score of 1 was assigned for each of the following conditions when they were true: Barthel score � 15 [58], MMSE score � 24 [60], FIM score � 100 [58] and residence setting in Community or Day-Hospital. Then, the subject was

classified as healthy if the sum of the scores was � 3 or unhealthy otherwise. Considering the elderly subjects from the ELDERMET study at T0 that were not outliers according to the total number of reads, we computed a predictive model for health status in old age based on GM biodiversity. After imputing missing data with the Multiple Imputation by Chained Equations (MICE) method [61] and standardizing the covariates, we computed a Principal Component Analysis (PCA) based on the elderly physical and cognitive state scores, (FIM score, Barthel index and MMSE score) and the residence setting (S5 Fig). Most of the variance (91.9%) was explained by the first Principal Component (PC[0]) that was hence con- sidered as the measure of healthy aging in the following predictive model. The predictive model of health status in old age based on GM biodiversity was computed as a Leave-One-Out cross validated linear regression and various possible GM biodiversity mea- surements were adopted as explanatory variables: besides θrare and θabundant, we also consid- ered the Pielou, Shannon, Simpson, Hill1, and Hill2 indices. After dichotomizing the subjects true and predicted health status in healthy/unhealthy, we computed the Area Under the Receiver Operating Characteristic Curve (AUC of ROC). The models were all adjusted for sex and standardized number of reads and computed both with and without the addition of calf circumference, BMI and the inflammatory markers IL-6, IL-8, IL-10, TNFa and CRP as furhter explanatory variables. Supporting information

S1 Fig. Model selection. Logarithm of the ratio of the posterior probabilities of model 2NB and 1NB (left), of model 2NB and 3NB (center) and of model 1NB and 3NB (right). Each pos- terior probability has been increased by a small factor (10−7) to avoid infinite value in the log- ratio. When comparing models 2NB and 1NB (left), data below the red line are those for which the 2NB model is selected rather than the 1NB. When comparing models 2NB and 3NB (center), data above the red line are those for which the 2NB model is selected rather than the 3NB. When comparing models 1NB and 3NB (right), data above the red line are those for which the 1NB model is selected rather than the 3NB. (EPS) PLOS ONE Gut microbiota biodiversity and healthy aging PLOS ONE | https://doi.org/10.1371/journal.pone.0237207 October 30, 2020 18 / 23

S2 Fig. Example of Preston plot and fit. Preston plot of the sample SRR3679961 of the Kong data set (gray histogram). The box-plots represent the predicted values obtained over the ABC iterations. The magenta and blue lines are two Negative Binomial distributions obtained using the median of the parameters of the ABC iterations respectively concerning “rare” and “abun- dant” species. (EPS) S3 Fig. AIC and BIC of natural cubic spline model. AIC and BIC (y-axis) of the natural cubic spline model adjsuted for sex and scaled total number of reads, when varying the num- ber of degrees of freedom (x-axis). Results for θrare are shown on the left and those for θabundant on the right. For both AIC and BIC, lower values are better. (EPS) S4 Fig. Biodiversity variability within subjects. For each subject of the ELDERMET data set for which data at the 3 time points were available, we plot the values of θrare and θabundant. Sub- jects are sorted according to the minimum average biodiversity. The plots show that the between-samples variability of both θrare and θabundant is higher than the within-sample vari- ability. (EPS) S5 Fig. Principal Component Analysis of the ELDERMET T0 health state. Representation of the ELDERMET elderly subjects (at T0) in the space defined by the first two components of the PCA computed using the FIM score, the Barthel index, the MMSE score and residence set- ting as covariates. Subjects are colored according to the healthy and unhealthy classification that was defined based on the same parameters used for the PCA. (EPS)

S6 Fig. Relationship between the variables of the ELDERMET T0 data set. Heat map of the Pearson’s r correlation matrix computed for the ELDERMET data set at T0. (EPS) S7 Fig. Results of the predictive model of health status in old age. Empirical and predicted values of the first PC. The plotted model was evaluated using as explanatory variables θrare and θabundant, plus calf circumference, BMI and the inflammatory markers IL-6, IL-8, IL-10, TNFa and CRP, and adjusting for sex and standardized total number of reads. AUC of ROC was computed after dichotomizing the subjects health status in healthy/unhealthy. Bisector repre- sents perfect prediction. AUC of ROC results for the alternative predictive models are reported in S2 Table. (EPS) S8 Fig. DS results obtained with classical biodiversity indices. (EPS) S9 Fig. ELDERMET results obtained with classical biodiversity indices. (EPS) S1 Table. Summary of the data sets. For each data set, we summarize the number of individu- als for each decade of age, divided by sex and health status. Here, we considered only subjects that were included in the statistical analysis. (XLSX) S2 Table. Summary of sequencing and bioinformatic processing. Title. For each data set, we report the main sequencing and bioinformatic methodologies applied to samples, including PLOS ONE Gut microbiota biodiversity and healthy aging PLOS ONE | https://doi.org/10.1371/journal.pone.0237207 October 30, 2020 19 / 23

(when available) sequencing platform and technology, reagent version, target region, algo- rithm used to merge reads (for paired reads), read filtering rules and OTU calculation. (XLSX) Author Contributions Conceptualization: Claudia Sala, Enrico Giampieri, Paolo Garagnani, Daniel Remondini, Armando Bazzani, Claudio Franceschi, Gastone C. Castellani. Data curation: Claudia Sala. Formal analysis: Claudia Sala, Enrico Giampieri. Funding acquisition: Daniel Remondini, Claudio Franceschi, Gastone C. Castellani. Methodology: Claudia Sala, Enrico Giampieri, Silvia Vitali, Armando Bazzani, Gastone C. Castellani. Software: Claudia Sala, Enrico Giampieri, Silvia Vitali. Visualization: Claudia Sala, Enrico Giampieri. Writing – original draft: Claudia Sala, Gastone C. Castellani. Writing – review & editing: Claudia Sala, Enrico Giampieri, Silvia Vitali, Paolo Garagnani, Daniel Remondini, Armando Bazzani, Claudio Franceschi, Gastone C. Castellani. References 1. Turnbaugh PJ, Hamady M, Yatsunenko T, Cantarel BL, Duncan A, Ley RE, et al. A core gut microbiome in obese and lean twins. nature. 2009 Jan; 457(7228):480–4. https://doi.org/10.1038/nature07540 PMID: 19043404 2. Marchesi JR, Adams DH, Fava F, Hermes GD, Hirschfield GM, Hold G, et al. The gut microbiota and host health: a new clinical frontier. Gut. 2016 Feb 1; 65(2):330–9. https://doi.org/10.1136/gutjnl-2015- 309990 PMID: 26338727 3. Claesson MJ, Jeffery IB, Conde S, Power SE, O’connor EM, Cusack S, et al. Gut microbiota composi-

tion correlates with diet and health in the elderly. Nature. 2012 Aug; 488(7410):178–84. https://doi.org/ 10.1038/nature11319 PMID: 22797518 4. Ferrucci L, Fabbri E. Inflammageing: chronic inflammation in ageing, cardiovascular disease, and frailty. Nature Reviews Cardiology. 2018 Sep; 15(9):505–22. https://doi.org/10.1038/s41569-018-0064-2 PMID: 30065258 5. Biragyn A, Ferrucci L. Gut dysbiosis: a potential link between increased cancer risk in ageing and inflammaging. The Lancet Oncology. 2018 Jun 1; 19(6):e295–304. https://doi.org/10.1016/S1470-2045 (18)30095-0 PMID: 29893261 6. Jeraldo P, Sipos M, Chia N, Brulc JM, Dhillon AS, Konkel ME, et al. Quantification of the relative roles of niche and neutral processes in structuring gastrointestinal microbiomes. Proceedings of the National Academy of Sciences. 2012 Jun 19; 109(25):9692–8. 7. Sala C, Vitali S, Giampieri E, do Valle I`F, Remondini D, Garagnani P, et al. Stochastic neutral modelling of the Gut Microbiota’s relative species abundance from next generation sequencing data. BMC bioin- formatics. 2016 Dec; 17(2):S16. https://doi.org/10.1186/s12859-015-0858-8 PMID: 26821617 8. Wilmanski T, Diener C, Rappaport N, Patwardhan S, Wiedrick J, Lapidus J, et al. Gut Microbiome Pat- tern Reflects Healthy Aging and Predicts Extended Survival in Humans. bioRxiv. 2020 Jan 1. 9. Lan Y, Kriete A, Rosen GL. Selecting age-related functional characteristics in the human gut micro- biome. Microbiome. 2013 Dec; 1(1):2. https://doi.org/10.1186/2049-2618-1-2 PMID: 24467949

10. Odamaki T, Kato K, Sugahara H, Hashikura N, Takahashi S, Xiao JZ, et al. Age-related changes in gut microbiota composition from newborn to centenarian: a cross-sectional study. BMC microbiology. 2016 Dec; 16(1):90. https://doi.org/10.1186/s12866-016-0708-5 PMID: 27220822 11. Galkin F, Aliper A, Putin E, Kuznetsov I, Gladyshev VN, Zhavoronkov A. Human microbiome aging clocks based on deep learning and tandem of permutation feature importance and accumulated local effects. bioRxiv. 2018 Jan 1:507780. PLOS ONE Gut microbiota biodiversity and healthy aging PLOS ONE | https://doi.org/10.1371/journal.pone.0237207 October 30, 2020 20 / 23

12. Mosca A, Leclerc M, Hugot JP. Gut microbiota diversity and human diseases: should we reintroduce key predators in our ecosystem?. Frontiers in microbiology. 2016 Mar 31; 7:455. https://doi.org/10. 3389/fmicb.2016.00455 PMID: 27065999 13. Xu C, Zhu H, Qiu P. Aging progression of human gut microbiota. BMC microbiology. 2019 Dec 1; 19 (1):236. https://doi.org/10.1186/s12866-019-1616-2 PMID: 31660868 14. Cardinale BJ, Duffy JE, Gonzalez A, Hooper DU, Perrings C, Venail P, et al. Biodiversity loss and its impact on humanity. Nature. 2012 Jun; 486(7401):59–67. https://doi.org/10.1038/nature11148 PMID: 22678280 15. Kreft JU, Plugge CM, Prats C, Leveau JH, Zhang W, Hellweger FL. From genes to ecosystems in micro- biology: modeling approaches and the importance of individuality. Frontiers in microbiology. 2017 Nov 27; 8:2299. https://doi.org/10.3389/fmicb.2017.02299 PMID: 29230200 16. Elmqvist T, Folke C, Nystro¨ m M, Peterson G, Bengtsson J, Walker B, et al. Response diversity, ecosys- tem change, and resilience. Frontiers in Ecology and the Environment. 2003 Nov; 1(9):488–94. 17. Loreau M, De Mazancourt C. Biodiversity and ecosystem stability: a synthesis of underlying mecha- nisms. Ecology letters. 2013 May; 16:106–15. https://doi.org/10.1111/ele.12073 PMID: 23346947 18. Larsen OF, Claassen E. The mechanistic link between health and gut microbiota diversity. Scientific reports. 2018 Feb 1; 8(1):1–5. https://doi.org/10.1038/s41598-017-17765-5 PMID: 29311619 19. Quercia S, Candela M, Giuliani C, Turroni S, Luiselli D, Rampelli S, et al. From lifetime to evolution:

timescales of human gut microbiota adaptation. Frontiers in microbiology. 2014 Nov 4; 5:587. https:// doi.org/10.3389/fmicb.2014.00587 PMID: 25408692 20. Biagi E, Franceschi C, Rampelli S, Severgnini M, Ostan R, Turroni S, et al. Gut microbiota and extreme longevity. Current Biology. 2016 Jun 6; 26(11):1480–5. https://doi.org/10.1016/j.cub.2016.04.016 PMID: 27185560 21. Cani PD, Delzenne NM. The role of the gut microbiota in energy metabolism and metabolic disease. Current pharmaceutical design. 2009 May 1; 15(13):1546–58. https://doi.org/10.2174/ 138161209788168164 PMID: 19442172 22. Franceschi C, Bonafè M, Valensin S, Olivieri F, De Luca M, Ottaviani E, et al. Inflamm-aging: an evolu- tionary perspective on immunosenescence. Annals of the new York Academy of Sciences. 2000 Jun; 908(1):244–54. https://doi.org/10.1111/j.1749-6632.2000.tb06651.x PMID: 10911963 23. Franceschi C, Garagnani P, Parini P, Giuliani C, Santoro A. Inflammaging: a new immune–metabolic viewpoint for age-related diseases. Nature Reviews Endocrinology. 2018 Oct; 14(10):576–90. https:// doi.org/10.1038/s41574-018-0059-4 PMID: 30046148 24. Furman D, Campisi J, Verdin E, Carrera-Bastos P, Targ S, Franceschi C, et al. Chronic inflammation in the etiology of disease across the life span. Nature medicine. 2019 Dec; 25(12):1822–32. https://doi. org/10.1038/s41591-019-0675-0 PMID: 31806905 25. Shannon CE. A mathematical theory of communication. Bell system technical journal. 1948 Jul; 27 (3):379–423. 26. Simpson EH. Measurement of diversity. Nature. 1949 Apr; 163(4148):688.

27. Pielou EC. The measurement of diversity in different types of biological collections. Journal of theoreti- cal biology. 1966 Dec 1; 13:131–44. 28. Hill MO. Diversity and evenness: a unifying notation and its consequences. Ecology. 1973 Mar; 54 (2):427–32. 29. Azaele S, Suweis S, Grilli J, Volkov I, Banavar JR, Maritan A. Statistical mechanics of ecological sys- tems: Neutral theory and beyond. Reviews of Modern Physics. 2016 Jul 26; 88(3):035003. 30. Hubbell SP. The unified neutral theory of biodiversity and biogeography (MPB-32). Princeton University Press; 2001 Apr 29. 31. Volkov I, Banavar JR, Hubbell SP, Maritan A. Patterns of relative species abundance in rainforests and coral reefs. Nature. 2007 Nov; 450(7166):45–9. https://doi.org/10.1038/nature06197 PMID: 17972874 32. Tang J, Zhou S. Hybrid niche-neutral models outperform an otherwise equivalent neutral model for fit- ting coral reef data. Journal of theoretical biology. 2013 Jan 21; 317:212–8. https://doi.org/10.1016/j.jtbi. 2012.10.019 PMID: 23089372 33. Biagi E, Candela M, Centanni M, Consolandi C, Rampelli S, Turroni S, et al. Gut microbiome in Down syndrome. PLoS one. 2014; 9(11). https://doi.org/10.1371/journal.pone.0112023 PMID: 25386941 34. Jeffery IB, Lynch DB, O’toole PW. Composition and temporal stability of the gut microbiota in older per- sons. The ISME journal. 2016 Jan; 10(1):170. https://doi.org/10.1038/ismej.2015.88 PMID: 26090993 35. Schnorr SL, Candela M, Rampelli S, Centanni M, Consolandi C, Basaglia G, et al. Gut microbiome of

the Hadza hunter-gatherers. Nature communications. 2014 Apr 15; 5:3654. https://doi.org/10.1038/ ncomms4654 PMID: 24736369 PLOS ONE Gut microbiota biodiversity and healthy aging PLOS ONE | https://doi.org/10.1371/journal.pone.0237207 October 30, 2020 21 / 23

36. Kong F, Hua Y, Zeng B, Ning R, Li Y, Zhao J. Gut microbiota signatures of longevity. Current Biology. 2016 Sep 26; 26(18):R832–3. https://doi.org/10.1016/j.cub.2016.08.015 PMID: 27676296 37. Bian G, Gloor GB, Gong A, Jia C, Zhang W, Hu J, et al. The gut microbiota of healthy aged Chinese is similar to that of the healthy young. Msphere. 2017 Oct 25; 2(5). https://doi.org/10.1128/mSphere. 00327-17 PMID: 28959739 38. Bazzani A, Sala C, Giampieri E, Castellani G. Master equation and relative species abundance distribu- tion for Lotka-Volterra models of interacting ecological communities. Theoretical biology forum.2016 Jan; 109(1–2):37–47. https://doi.org/10.19272/201611402003 PMID: 29513351 39. An R, Wilms E, Masclee AA, Smidt H, Zoetendal EG, Jonkers D. Age-dependent changes in GI physiol- ogy and microbiota: time to reconsider?. Gut. 2018 Dec 1; 67(12):2213–22. https://doi.org/10.1136/ gutjnl-2017-315542 PMID: 30194220 40. Garcı´a-Peña C, A´ lvarez-Cisneros T, Quiroz-Baez R, Friedland RP. Microbiota and aging. A review and commentary. Archives of medical research. 2017 Nov 1; 48(8):681–9. https://doi.org/10.1016/j.arcmed. 2017.11.005 PMID: 29229199 41. Jackson MA, Jeffery IB, Beaumont M, Bell JT, Clark AG, Ley RE, et al. Signatures of early frailty in the gut microbiota. Genome Med. 2016 Jan 29; 8(1):8. https://doi.org/10.1186/s13073-016-0262-7 Erratum in: Genome Med. 2016; 8(1):21. Jackson, Matt [corrected to Jackson, Matthew A]. PMID: 26822992; PMCID: PMC4731918.

42. Biagi E, Nylund L, Candela M, Ostan R, Bucci L, Pini E, et al. Through ageing, and beyond: gut micro- biota and inflammatory status in seniors and centenarians. PloS one. 2010; 5(5). 43. Maffei VJ, Kim S, Blanchard IV E, Luo M, Jazwinski SM, Taylor CM, et al. Biological aging and the human gut microbiota. Journals of Gerontology Series A: Biomedical Sciences and Medical Sciences. 2017 Oct 12; 72(11):1474–82. https://doi.org/10.1093/gerona/glx042 PMID: 28444190 44. Gordon A, Hannon GJ. Fastx-toolkit. FASTQ/A short-reads preprocessing tools (unpublished) http:// hannonlab.cshl.edu/fastx_toolkit. 2010 Jan 20; 5. 45. Edgar RC. UPARSE: highly accurate OTU sequences from microbial amplicon reads. Nature methods. 2013 Oct; 10(10):996. https://doi.org/10.1038/nmeth.2604 PMID: 23955772 46. Schloss PD, Westcott SL, Ryabin T, Hall JR, Hartmann M, Hollister EB, et al. Introducing mothur: open- source, platform-independent, community-supported software for describing and comparing microbial communities. Appl. Environ. Microbiol. 2009 Dec 1; 75(23):7537–41. https://doi.org/10.1128/AEM. 01541-09 PMID: 19801464 47. Raftery AE. Bayes factors and BIC: Comment on “A critique of the Bayesian information criterion for model selection”. Sociological Methods & Research. 1999 Feb; 27(3):411–27. 48. Marin JM, Pudlo P, Robert CP, Ryder RJ. Approximate Bayesian computational methods. Statistics and Computing. 2012 Nov 1; 22(6):1167–80. 49. Python Software Foundation. Python Language Reference, version 3.6. Available at http://www.python. org

50. patsy—Describing statistical models in Python. Available at https://patsy.readthedocs.io/en/latest/ 51. McKinney W. Data structures for statistical computing in python. In Proceedings of the 9th Python in Science Conference. 2010 Jun 28;445:51–56. 52. Walt SV, Colbert SC, Varoquaux G. The NumPy array: a structure for efficient numerical computation. Computing in Science & Engineering. 2011 Mar; 13(2):22–30. 53. Virtanen P, Gommers R, Oliphant TE, Haberland M, Reddy T, Cournapeau D, et al. SciPy 1.0: funda- mental algorithms for scientific computing in Python. Nature methods. 2020 Mar; 17(3):261–72. https:// doi.org/10.1038/s41592-019-0686-2 PMID: 32015543 54. Hunter JD. Matplotlib: A 2D graphics environment. Computing in science & engineering. 2007 May; 9 (3):90–5. 55. Waskom M, Botvinnik O, Hobson P, Warmenhoven J, Cole JB, Halchenko Y, et al. Seaborn: statistical data visualization. URL: https://seaborn.pydata.org/. 2014. 56. Patil A, Huard D, Fonnesbeck CJ. PyMC: Bayesian stochastic modelling in Python. Journal of statistical software. 2010 Jul; 35(4):1. PMID: 21603108 57. Hastie TJ, Tibshirani RJ. Generalized additive models. CRC press; 1990 Jun 1. 58. Granger CV, Hamilton BB, Keith RA, Zielezny M, Sherwin FS. Advances in functional assessment for medical rehabilitation. Topics in geriatric rehabilitation. 1986 Apr 1; 1(3):59–74. 59. Mahoney FI, Barthel DW. Functional evaluation: the Barthel Index: a simple index of independence use- ful in scoring improvement in the rehabilitation of the chronically ill. Maryland state medical journal. 1965. PLOS ONE

Gut microbiota biodiversity and healthy aging PLOS ONE | https://doi.org/10.1371/journal.pone.0237207 October 30, 2020 22 / 23

60. Pangman VC, Sloan J, Guse L. An examination of psychometric properties of the mini-mental state examination and the standardized mini-mental state examination: implications for clinical practice. Applied Nursing Research. 2000 Nov 1; 13(4):209–13. https://doi.org/10.1053/apnr.2000.9231 PMID: 11078787 61. White IR, Royston P, Wood AM. Multiple imputation using chained equations: issues and guidance for practice. Statistics in medicine. 2011 Feb 20; 30(4):377–99. https://doi.org/10.1002/sim.4067 PMID: 21225900 PLOS ONE Gut microbiota biodiversity and healthy aging PLOS ONE | https://doi.org/10.1371/journal.pone.0237207 October 30, 2020 23 / 23