Global and temporal state of the human gut microbiome in health and disease

The role of gut microbiota in humans is of great interest, and metagenomics provided the possibilities for extensively analysing bacterial diversity in health and disease. Here we explored the human gut microbiome samples across 19 countries, performing compositional, functional and integrative analysis. To complement these data and analyse the stability of the microbiome, we followed 86 healthy Swedish individuals over one year, with four sampling times and extensive clinical phenotyping. The integrative analysis of temporal microbiome changes shows the existence of two types of species with a tendency to vary in abundance with time, here called outow and inow species. Importantly, the former tends to be enriched in disease, while the latter is enriched in health. We suggest that the decrease of disease-associated outow and the increase of health-associated inow species with time may be a fundamental albeit previously unrecognized aspect of the homeostasis maintenance in a healthy microbiome. researchers to explore for the rst time an integrative analysis on composition, functional, richness, diseases and region signatures for the gut microbiota across 19 geographical regions and 20 diseases. This analysis was followed by investigating the gut microbiome of healthy Swedish individuals with four times sampling across one year to study the longitudinal variability of the microbiota. This revealed the tendency of disease-associated species to decrease in abundance. In contrast, the health-associated species tended to increase in abundance. We suggest that this dynamic contributes to the maintenance of gut microbial homeostasis in healthy individuals. These ndings were further validated by follow up sampling from same cohort with additional two time points across 6 months. of three based on F-tests of regressions (p-value < 0.01; β Ei , β Ni , β Ci > 0) and quaternary plots were shown based on squared regression coecients (β Ei2 , β Ni2 , β Ci2 ) normalized by their total sum. one-sided in countries such as A. histaminiformans (NAFLD), and species depleted, such as G. urolithinfaciens (NAFLD). h, Jitter plots of frequency of the signicantly enriched/depleted cohorts of all MSPs were calculated: total frequency of enriched/depleted cohorts (|number of enriched cohorts| + |number of depleted cohorts|) and subtracted frequency between enriched cohorts and depleted cohorts (|number of cohorts| of depleted Point colours Note:


Introduction
Metagenomic studies of the human microbiome enable the characterization of the microbial and functional diversity in health and disease 1 . Advances in metagenome assembly and various clustering methods enabled the generation of metagenome species [2][3][4][5][6] . Most of these studies focused on unveiling new uncultured genomes, while only few focused on investigating the functional potentials and dynamic changes of the gut microbiome [7][8][9] . Understanding the functional and temporal behaviour of the microbiome may have great implications for the identi cation of its global signature in health and disease [10][11][12] . Additionally, short-term perturbations may trigger gut microbiota dysbiosis and changes at compositional and functional levels. Speci cally, the negative selective microbe-microbe and hostmicrobe interactions, in the context of metabolism or antimicrobial machinery, could be the main mechanism underlying microbial dysbiosis 13 . Large-scale integration of microbiome functional changes and their associations with clinical data may provide novel information on temporal changes in the microbiome and host physiology 14 .
Herein, we integrated publicly available data from a large number of studies across different countries from healthy and diseased individuals. The analysis is presented in an open-access Human Gut Microbiome Atlas (www.microbiomeatlas.org), allowing researchers to explore for the rst time an integrative analysis on composition, functional, richness, diseases and region signatures for the gut microbiota across 19 geographical regions and 20 diseases. This analysis was followed by investigating the gut microbiome of healthy Swedish individuals with four times sampling across one year to study the longitudinal variability of the microbiota. This revealed the tendency of disease-associated species to decrease in abundance. In contrast, the health-associated species tended to increase in abundance. We suggest that this dynamic contributes to the maintenance of gut microbial homeostasis in healthy individuals. These ndings were further validated by follow up sampling from same cohort with additional two time points across 6 months. Human Gut Microbiome Atlas; Pan-metagenomics study on compositional and functional changes of human gut microbiome We performed large-scale integrative analysis of 4,880 publicly available shotgun metagenomics stool samples, with at least 10 million high-quality sequencing reads from healthy and diseased cohorts from 19 different countries across ve continents ( Fig. 1a-b and Supplementary Table S1). We rare ed all metagenomic samples into 10 million reads per sample, which enable comparative analysis across different cohorts. We created the Human Gut Microbiome Atlas (HGMA) using quantitative analysis of shotgun metagenomics based on microbial genomes assembled using Metagenomic Species Pangenomes (MSPs) (Fig. 1c). The MSP number was increased from 1,661 (previous release 5 ) to 1,989 (average number of genes 1,894 ± 1,616) (Methods), and their taxonomy was updated. We generated gene counts and MSP abundances for all the samples using the 10.4 million gene catalogue 15 . We also characterized the functions and phenotype of the MSPs in 7 different categories (KO, PFAM, CAZyme, Mustard, JGI-GOLD phenotype, PATRIC virulence factor, and antiSMASH biosynthetic gene clusters) and identi ed co-conserved functional clusters across species (7,763 clusters) (Methods). This information was completed with 344 newly sequenced longitudinal samples from 86 Swedish individuals, described in detail in a subsequent section (Fig. 1d). All the data are freely available in the HGMA, without restrictions, in the public open access database (www.microbiomeatlas.org) that is part of the Human Protein Atlas program (https://www.proteinatlas.org). All MSPs and functions are highlighted together with the 5,224 samples across 19 countries with disease and healthy cohorts.
Using the 3,039 samples obtained from healthy individuals across 18 countries, including westernized and non-westernized regions, we uncovered the geographical distribution of the healthy gut microbiome. To this end, we applied the unsupervised clustering method, monocle, to MSP abundance pro les of the 3,039 samples (Methods) 16,17 . We observed that there were two distinct ordinations of non-westernized and European samples of healthy subjects (  Table S2). Functional annotation-based analysis across geographical clusters, revealed an enrichment of CAZymes for degrading host mucins and storage carbohydrates in westernized population, where antimicrobial resistance (AMR) and virulence factors were also more prevalent (Fig. 1f). Comparison of functions of region-enriched MSPs in European countries and in US/China/Japan revealed that fosfomycin and aminoglycoside resistance were the signi cant AMR, as deduced from the explained variances of ANOVA (Methods). Among the biosynthetic genes encoding the production of secondary metabolites, resorcinol, lanthipeptide, bacteriocin and homoserine lactone were enriched in the European and US/China/Japan populations. These secondary metabolites, together with AMR, appeared to be the key feature in the westernized countries.
To distinguish diseased and healthy microbiome from multiple cohorts, we performed a panmetagenomics association study (Pan-MGAS) of multiple disease cohorts (18 diseases across 28 cohorts from 11 westernized countries). We determined the enriched and depleted species in diseased compared to healthy control samples, with an effect size > 0.3 ( Fig. 1g and Supplementary Table S3). Speci c species were either enriched or depleted in a certain disease(s), regardless of geographical differences. For instance, in individuals with fatty liver disease, Gordonibacter urolithinfaciens and Allisonella histaminifomans were depleted and enriched, respectively. We also found that species associated with low gene richness 18 , such as Clostridium bolteae, were enriched in several diseases. Similarly, Parvimonas micra was enriched in colorectal cancers, regardless of geographical differences. Furthermore, Akkermansia muciniphila was often depleted in several diseases. A. muciniphila is associated with improved intestinal barrier function and its depletion suggests intestinal barrier disruption in these different diseases 19 . The analysis of the frequency of the enriched/depleted MSPs among diseased cohorts showed that there were common species that could initially disrupt the microbiome balance and cause gut dysbiosis or foster microbial symbiosis that promotes health (top-left and top-right boxes in Fig. 1h, Supplementary Table S4).
Interestingly, we observed that MSPs commonly depleted in diseases were highly country-speci c, while commonly enriched MSPs were usually shared by several diseases and were less related to geographical origins (Extended Fig. 1e-i and Supplementary Table S4). Moreover, we observed that MSPs commonly enriched or depleted in diseases were associated with different temporal behaviours among healthy individuals, as detailed in the following sections.
Dynamic changes of gut microbiome composition; in ow and out ow species We investigated the temporal changes in microbiome composition at the individual level by applying Markov chain models (MCMs) to the MSPs identi ed in the longitudinal cohort of 86 healthy Swedes (Methods). This analysis enabled us to estimate the transition probability of individual MSPs from presence to absence (out ow probability) and vice versa (in ow probability) across different sampling times. We identi ed two groups of MSPs preferably transiting from presence to absence or from absence to presence; for brevity, we term them "out ow species (OFS)" and "in ow species (IFS)" respectively (  Table S5). Clearly, declaring a species absent or present depends on the detection threshold, which is in turn determined by sequencing depth. We performed the analysis at three depth levels for 15, 10, and 5 million reads, and observed largely concordant results (Extended Fig. 3). For instance, 35 IFS (90%) were detected at both 10 and 15 million reads levels, while 4 and 6 species were detected only at former and latter, respectively. Similar results were observed for OFS: 447 (88%) were detected at both levels, while 62 and 27 species were detected at 10 and 15 million reads only.
Overall, in ow and out ow probabilities were highly correlated at three different depth levels, with a slight reduction for out ows at 5 million reads (Supplementary Table S5).
To determine the robustness of these ndings, we compared species-retaining probability (Kaplan-Meier estimates) with out ow probability, expected to be inversely proportional (Extended Fig. 4a and Methods). We observed that the species that had lower out ow probability, such as Bacteroides vulgatus and Prevotella copri, indeed had the highest retaining probability, whereas those of higher out ow probability (e.g. Veillonella infantium and Streptococcus parasanguinis) had reduced retaining probability. The species retaining probabilities were correlated with their mean abundance (Extended Fig. 4b), even if the associations did not appear signi cant for any individual species based on Cox regression (p-values > 0.1, Extended Fig. 4c).
We observed that the changes of OFS abundances among 86 individuals (Δ OFS ) were inversely correlated with those of IFS (Δ IFS ), i.e. suggesting competition between IFS and OFS (Spearman's correlation = -0.334, p-value = 4.6 × 10 -8 ; Fig. 2b-c, Extended Fig. 2g). A higher abundance of OFS increased the gut microbiome imbalance between different visits, as the similarity of the gut microbiome compositions between the visits was decreased in 31 OFS-enriched individuals compared to 37 OFS-depleted individuals (Wilcoxon one-sided test p-value < 0.01; Fig. 2d). IFS-enriched individuals maintained their gut microbial composition such that it was similar between different time points (Extended Fig. 2h). Interestingly, increasing abundance of IFS was associated with increasing gene richness, known to be related to better health 18 , suggesting that IFS may be bene cial (Spearman's correlation = 0.206, p-value = 9.0 × 10 -4 ; Fig. 2e and Extended Fig. 2i).
We hypothesised that IFS and OFS may differ in their growth rates, the former outgrowing the latter, and tested this hypothesis in three ways. First, we estimated species growth rates from metagenomic samples by Growth Rate InDex (GRiD) analysis 20 (Methods). For that, we strati ed individuals into groups, enriched in IFS or OFS and found that in both groups GRiD scores of IFS were higher than OFS; they were the highest among IFS-enriched group ( Fig. 2f-g). Second, we assessed species growth rates in bioreactors inoculated with healthy human stool samples, via GRiD analysis ( Fig. 2h-i, Methods). We observed that the growth of IFS increased signi cantly over 24 hours, whereas that of OFS did not change, demonstrating that IFS could outgrow the OFS. Third, we used genome scale modelling to simulate species growth rates 21-24 (Methods) on four different putative diets (high bre and high protein for plant-and animal-based diets) for highly prevalent 34 OFS and 30 IFS. The predicted growth rates of the IFS were signi cantly higher than of OFS ( Fig. 2j-k, Supplementary Table S6). Furthermore, the reaction essentiality analysis indicated that the out ow GEMs were signi cantly dependent on the substrate, often displaying amino acid auxotrophy (Supplementary Figure 3 and 4). We suggest that the differences in growth rates and substrate dependence between IFS and OFS could underlie the directionality of the gut microbiome dynamics we report.
To test whether the in ow and out ow assignments of species deduced from the analysis of the four time points in our cohort persist over time, we collected and analysed two additional time points with 6 months intervals from the same individuals ( Fig. 2l-m). Furthermore, to examine whether the assignments de ned from a Swedish study are also found in other, geographically different regions, we analysed two publicly available cohorts, from Italy and US 9,25 ( Fig. 2n-q). In all cases, for both IFS and OFS, the transition probabilities were well correlated with those found for the rst time points of our longitudinal cohort (Spearman's correlation coe cients > 0.56) for all comparisons). We conclude that IFS and OFS are largely conserved and are thus a global feature of the human gut microbiome.
Enrichment of IFS and OFS species determines richness, dysbiosis, and host physiology To further explore links between gene richness and IFS and OFS, we determined the gene richness distribution of HGMA samples and grouped healthy samples as either high or low gene richness (HGR and LGR) based on the top and bottom 25% gene counts of HGMA samples, respectively (Extended Fig.  5a-c). We then identi ed species enriched in HGR and LGR by comparing MSP abundances: total 759 MSPs were enriched in HGR and 95 MSPs were enriched in LGR (|Log 2 fold change| > 2, Wilcoxon rank two-sided test adjusted p-value < 10 -3 ) (Supplementary Table S7). Interestingly, LGR-enriched MSPs were signi cantly enriched in OFS (i.e. higher out ow probability), while no enrichment was observed for HGRenriched MSPs (Fig. 3a-b). Low gene richness was previously associated with a risk of developing chronic diseases related to the metabolic syndrome; the enrichment of disease associated OFS in LGR individuals was thus coherent with that observation 18 . We observed similar species enriched with low gene richness and associations with metabolic phenotypes (Extended Fig. 5d).
To investigate the impacts of IFS and OFS species on host physiology, we next examined the association of IFS and OFS with clinical parameters among healthy individuals (Supplementary Table S8). We traced the changes in 40 clinical parameters of 86 Swedish individuals (Supplementary Table S9) and linked them with the abundances of IFS and OFS species using linear mixed effect models (Methods). We found that IFS abundances were associated with muscle composition (p-value = 0.018, Fig. 3c), thus showing associations with microbial core metabolism such as amino acid metabolism, whereas OFS abundances were associated with a more diverse spectrum of clinical parameters, such as blood glucose, urate, B-type natriuretic peptide (BNP), ApoA1, and erythrocyte particles (p-value < 0.05). Interestingly, the OFS abundances were positively associated with BNP as a key heart failure marker (Extended Fig. 6c-g, Supplementary Table 10).
Next, we associated the common depleted and enriched species in diseases (Fig. 1h), from Pan-MGAS analysis of 18 diseases, to the IFS and OFS. The commonly depleted MSPs had a greater tendency to be IFS ( Fig. 4a-b). On the contrary, the commonly enriched species (top-right, Fig. 1h) were likely to be OFS, compared to the depleted species across all diseases. Out of the 23 commonly enriched species in diseases (enriched in at least 3 different disease cohorts), 14 (61%) were OFS species, signi cantly more than the OFSs species within all species detected in the cohorts (36%) (Chi-square test, p-value = 0.015, Supplementary Table S4). In addition, among the 14 OFS species, 5 (36%) were opportunistic pathogens reported previously 26 (Chi-square test, p-value = 10 -9 ). These observations support the view of microbiome dynamics lowering the abundance of potentially harmful species in healthy adults.

Functional understanding of region-enriched species, IFS and OFS
Our functional analysis indicated that the IFS species were enriched in core metabolism, essential for energy homeostasis and for biosynthesis of macromolecules (i.e., amino acids, carbohydrates and fatty acids; Fig. 4c and Supplementary Table S11 and S12, Methods). They were also enriched in: (i) processes associated with increased survival, such as sporulation, cobalamin biosynthesis (CobS), and sirohydrochlorin cobaltochelatase (CbiK); (ii) secondary metabolites (bacteriocin); (iii) proteins related to starch and plant-based bre use (CAZymes GT5, GH13, GH51); (iv) anaerobic phenotypes (Supplementary Table S11). By contrast, the OFS species were enriched in accessory metabolism, such as biodegradation of xenobiotics (benzene, toluene, ethylbenzene, and xylenes -BTEX), paralleled by that of ABC transporters, possibly involved in the import of xenobiotics, suggesting that exposure to pollutants may promote their appearance (Fig. 4d, Extended Fig. 6a-b, and Supplementary Table S11). They were also enriched in (i) active sugar transport (i.e., phosphotransferase system (PTS); (ii) virulence factors (VFs) and trigger factors; (iii) putative competence protein ComGF and type IV secretion systems, the latter two being important mechanisms for horizontal gene transfer 27 . Finally, those microbes tended to be facultative anaerobes and microaerophiles.
Out ow enriched-functional clusters showed distinct links to gut microbiome dysbiosis To better understand potential impact of IFS and OFS at the functional module level, we identi ed co-conserved functional clusters of microbiome by applying an unsupervised clustering approach on MSPs (Fig. 4e, Extended Fig. 7 and Methods). This analysis provided a better representation of microbial functions than single annotations or known pathway de nitions (e.g. KEGG) (Extended Fig. 8). From the community detection algorithm, we identi ed 7,763 functional clusters, 6,297 singletons, and 591 representative clusters (Methods, Supplementary Table 13). For example, antimicrobial resistance and secondary biosynthetic genes were found to be singletons and not co-conserved with other functional genes. After excluding singletons and unreliable functional clusters detected in less than three species, we retained 591 representative clusters of microbial functions. One of the two largest clusters (CL-12 in Supplementary Table 13, named "comm-cluster" herewith) was over-represented among many commensal species, while the other (CL-10, named "patho-cluster") was enriched in a few pathobionts, such as Klebsiella spp., Enterobacter spp., and E. coli. Interestingly, the comm-cluster was enriched with genes involved in the biosynthesis of amino acids indicative of functions enriched in IFS species. In contrast, the patho-cluster was enriched in functions associated with the uptake of several substrates, again indicative of transporters enriched in OFS species. These included siderophore, ion, amino acid, and vitamin transport, thus competing with host and commensal bacteria. We also found other enrichedfunctional clusters, such as butyrate metabolism, propionate metabolism, vitamin B12, coenzyme metabolism, chemotaxis, ATPase, and mobile genetic elements (i.e., integrase and transposase) and the CRISPR-cas system (Fig. 4e); a number of these were correlated with phylum-level taxonomy (Extended Fig. 7c).
We next projected the functional clusters on enriched/depleted MSPs in HGMA disease cohorts ( Fig. 4f and Supplementary Fig. 5: hypergeometric tests, p-value < 10 -3 ). We found that many of disease-enriched functional clusters were enriched in the OFS species, for example, isoprenoid biosynthesis, competence proteins for DNA transformation, key signatures of OFS species, virulence factors, and nutrient uptake (e.g. ascorbate and mannose). It has been previously shown that isoprenoid biosynthesis initiates the majority of secondary metabolism 28 . However, we found a few functional clusters associated with species depleted in diseases, such as the CRISPR-cas system (i.e., the bacterial immune system) and teichoic acid transport.

Discussion
We have performed a comprehensive integrative analysis of global and temporal gut microbiomes, and we provide an open access HMGA portal (http://microbiomeatlas.org). Con rming previous observations 29 , we have described the gut microbiome regional speci city, which needs to be taken into account before using the gut microbiome for the strati cation of patients or for designing intervention studies. Beyond previous observations, our function-based analysis indicates that the western-enriched bacteria might dominate the gut microbial community with the production of antimicrobial peptides and homoserine lactone, which may inhibit their competitors.
Previous studies reported the temporal stability of the gut microbiome composition in an individual 7-9 , implying oscillations around an average value. Our integrative analysis of temporal microbiome changes in a longitudinal study of healthy individuals has shown the existence of directionality of compositional variations: there are two types of species with a tendency to either increase or decrease in abundance with time, termed in ow or out ow species, respectively. Importantly, out ow species include most of the known opportunistic pathogens, while in ow species are essentially devoid of them. Remarkably, our function-based analysis indicates that out ow species might have a negative impact on host physiology, as they have enriched accessory metabolism and secretion of virulence factors. Most interestingly, out ow species tend to be enriched in different diseases while, in contrast, the in ow ones tend to be enriched in healthy individuals. We suggest that the tendency for the former to decrease and the latter to increase in healthy individuals is a previously unrecognized facet of the gut microbiome homeostasis.
The out ow species tend to be facultative anaerobes and to have an oral origin (e.g., Streptococcus spp.). This observation suggests an increase in oral microbial transmission to the gut, possibly due to a decrease in gut microbiome resilience. Enrichment of oral species in the gut has been observed in several diseases 30,31 , and we suggest that increased mouth to gut microbial ow could be one of the global features of dysbiosis.
We have described the temporal dynamics of the gut microbiome through the discovery of out ow and in ow species. The enrichment of in ow species in healthy populations could possibly be due to their involvement in the storage carbohydrate degradation, such as starch and bre, accounting for their higher persistence. We consider two mechanisms of out ow species enrichment in disease. First, the out ow species were enriched in competence mechanisms, facilitating import of genetic elements such as AMR, possibly conferring selective advantage in the gut. Enrichment of drug e ux mechanisms in out ow species might also confer resistance to antibiotics and other medications used in disease treatments. Second, the out ow species may plunder nutrient uptake from the host and utilize simple monosaccharides to increase their abundance. Metabolic modelling indicates that IFS are favoured by common diets. It is known that short-chain fatty acids (SCFAs) have an important impact on microbemicrobe interactions and host physiology. However, we also observed that the comm-cluster that comprises the biosynthesis of amino acids and folate metabolism may have a larger impact on the resilience of the gut microbiome and host physiology. We observed that the enrichment of the patho-cluster could gradually or instantly perturb the gut ecosystem and shift it to either short-term gut imbalance or a persistent dysbiotic state.
Finally, the integration of metagenomics data from a large number of studies spanning ve continents provides valuable knowledge for researchers interested in the impact of the microbiome on individual health parameters. The open-access atlas will be updated routinely with the new publicly available gut metagenomics data, including the recently announced one million microbiome project aimed at providing comprehensive open-access metagenomics data from multiple research centres. In this manner, in-depth analysis of the impact of the gut microbiome on health and disease will be used to facilitate future studies to reveal the key role of the gut microbiome in human maintaining health.

Methods
Metagenomics species pan-genome (MSP) creation 1601 metagenomic samples used to build the Integrated Gene Catalog of the human gut microbiome (IGC2) were downloaded from the European Nucleotide Archive. Using the Meteor software suite 1 , reads from each sample were mapped against the IGC2 catalog and a raw gene abundance table was generated. This table was submitted to MSPminer 2 that reconstituted 1,989 clusters of co-abundant genes named Metagenomic-Species Pangenomes (MSPs). Quality control of each MSP was manually performed by visualizing heatmaps representative of the normalized gene abundance pro les. In addition, MSPs completeness and contamination was assessed by searching for 40 universal single copy marker genes 3 and by checking taxonomic homogeneity.

MSP taxonomic annotation with phylogenetic tree
MSPs taxonomic annotation was performed by aligning all core and accessory genes against nt and NCBI WGS (version of September 2018 restricted to the taxa Bacteria, Archaea, Fungi, Viruses and Blastocystis) using blastn (version 2.7.1, task = megablast, word_size = 16) 4 . The 20 best hits for each gene were kept. A species-level assignment was given if more than 50% of the genes matched the RefSeq reference genome of a given species, with a mean identity ≥ 95% and mean gene length coverage ≥ 90%. The remaining MSPs were assigned to a higher taxonomic level (genus to superkingdom), if more than 50% of their genes had the same annotation. 40 universal phylogenetic markers genes were extracted from the MSPs with 5 . MSPs with less than 5 markers were discarded. Then, the markers were separately aligned with MUSCLE 6 . The 40 alignments were merged and trimmed with trimAl 7 . Finally, the phylogenetic tree was computed with FastTreeMP 8 and visualized with iTOL 9 .
Phylogenetic placement was further used to improve and correct taxonomic annotation.
Wellness study population, sample collection, extraction, library prep and sequencing The wellness study is an ongoing prospective cohort study based on the Swedish CArdioPulmonary bioImage Study (SCAPIS) with 30,154 individuals enrolled at ages between 50 and 64 years recruited from random sampling of the general Swedish population. A total of 101 healthy individuals were recruited in the study and followed longitudinally for two years. Examinations in SCAPIS include imaging to assess coronary and carotid atherosclerosis, clinical chemistry, anthropometry, and extensive questionnaires, as previously described 10 . All participants provided written informed consent. The study protocol conforms to the ethical guidelines of the 1975 Declaration of Helsinki.
Total genomic DNA was isolated from 100-120 mg of faces using a repeated bead beating method.  11 . Protein sequences were aligned against 9,462 ARD sequences using blastp 2.7.1+ (option -evalue = 10 -5 ). Best-hit alignments were ltered for identity ≥ 95% and bidirectional alignment coverage ≥ 90% (at query and subject level), giving a list of ARD candidates belonging to 30 families. Annotation of the carbohydrate-active enzymes (CAZymes) of the IGC2 catalog was performed by comparing the predicted protein sequences to those in the CAZy database and to Hidden Markov Models (HMMs) built from each CAZy family 12 , following a procedure previously described for other metagenomics analysis 13 . Proteins of IGC2 catalog were also annotated to KEGG orthologous using Diamond (version 0.9.22.123) 14 against KEGG database (version 82). Best-hit alignments with e-value ≤ 10 -5 and bit score ≥ 60 were considered. Proteins involved in virulence factors of PATRIC 15,16 were matched against IGC2 17 by BLASTP (best identity > 50%, e-value < 10 -10 ). Phenotype of MSP were manually checked and annotated based on JGI-GOLD phenotype (organism metadata) 18 .
We identi ed biosynthetic genes of MSP with the use of standalone anti-SMASH program with minimal run option, focused on core detection modules (version 5) 19 . Loading antiSMASH into Amazon cloud computing (AWS) as docker image, we executed its mining process per MSP in a massive parallel setting.

Quality control/normalization of gene counts and species abundance pro ling
We ltered out human reads and then mapped metagenomic data (Supplementary Table 1) on IGC2 catalogue of human gut metagenome by METEOR 1 and based on the aligned reads, we estimated the abundance of each reference gene of the catalogue, normalizing multiple mapped reads by their numbers and summing up normalized counts for a given gene. Reducing the variability by sequencing depths, gene count values were downsized into 10 million reads per sample; and any samples less than 10 million mapped reads were excluded from our dataset. Normalized gene counts were used for the quanti cation of MSP abundance by R momr (MetaOMineR) package 20,21 . MSP abundances were estimated by the median abundance of the 25 marker genes representing the robust centroid of gene clusters of MSP. Sample metadata of all metagenomics data such as sequencing platform, geography, age, body-mass index, gender and the data source were provided under HGMA (http://microbiomeatlas.org).

Tracing the diversi cations of healthy metagenomic samples of different geography
After the quanti cation and per-million scaling of MSP abundance pro les, we employed trajectory analysis in R monocle ver.2 package to identify how samples were clustered 22 . In short, we selected the species pro les of all normal samples from different geographical origins and reduced the sample pro les into two dimensions by advanced nonlinear reconstruction algorithm, DDRTree. Based on the reduced two-dimensional components, we presented how samples were closely clustered as branches in scatter plots. Based on reduced pro les, we also calculated centroids and standard deviations of samples of given countries, except Finland population in toddlers (2 years).

Identi cation of region-enriched species from geographically distinct cohorts
We selected healthy samples of 17 countries after excluding matched controls of two-year old subjects of Finland T1D cohort and redundant samples of subjects with multiple measurements (i.e. multiple visits).
Among 17 countries, we estimated effect sizes for Wilcoxon signed rank (one-sided) tests 23 of different MSP abundances of two different countries. As one-sided tests were used, we set the lower bound of effect sizes as zero and the upper bound of effect sizes as one, avoiding negative and in nite values.
Based on estimated effect sizes, we identi ed signi cantly enriched species having medium effect sizes of speci c country (effect size ≥ 0.3), compared to six or more countries, and de ned those species as "region-enriched" species.

Next we categories species if enriched in 1) European countries, 2) non-westernized countries, and 3)
China/Japan/US and identi ed contrasted functions among those three clusters of countries by multivariate regressions as follows: where Y i indicates a function regard to species i like CAZyme, antibiotics resistance, anti-SMASH, and virulence factor (if a given function exists in species i, Y i =1, otherwise Y i = 0), ϵ indicates an intercept, For the validation of in ow and out ow from same Swedish wellness cohort, we additionally followed the two more visits (by every three months) and processed metagenomics data of 67 subjects (134 samples) after excluding subjects of missing visits and low sequencing depth less than 10 million reads. For the validation of in ow and out ow from independent cohorts, we processed metagenomics data from Italy (DINAMIC cohort) and US (HPFS cohort) after excluding subjects of missing visits and low sequencing depth less than 10 million reads. In HPFS cohort, we only took six-months interval samples of individuals, excluding one-day interval samples. We counted presence/absence of MSPs from the abundance pro les in a similar way of calculation in Swedish wellness cohort, and calculated state transition probabilities between presence and absence (i.e. in ow and out ow) after tting presence/absence pro les into twostate Markov chain model. Associations between MSP abundance pro les and clinical metadata Scaled abundance of IFS and OFS species populations together (Z IFS and Z OFS , respectively) were tested their associations with clinical parameters with considering random effects of individuals by linear mixed-effect models using R lme4 packages (p-values < 0.05) like below: where Y is clinical parameter, β IFS and β OFS are coe cients of xed effect variables, Z IFS and Z OFS , respectively, u i is a random intercept for subject i, and ϵ is residual.
In addition, we tested associations of single MSP with clinical parameters of given samples of wellness cohorts by linear mixed effect models like below: where Y is clinical parameter, β i is coe cient of xed effect variable, A i , u j is a random intercept for subject j, and ϵ is residual. We identi ed signi cant associations between MSP abundance and clinical parameters based on explained variance of xed effect calculated using R MuMIn package (explained variance > 10%).

Fecal fermentation in ARCOL bioreactor
M-ARCOL is a one-stage fermentation system run under semi-continuous conditions that simulates the main physicochemical and microbial conditions encountered in the human colonic ecosystem 25 . It consists of pH and temperature controlled, stirred (400 rpm), airtight glass vessels inoculated with fecal samples from human volunteers and maintained under anaerobic conditions by the sole activity of resident microbiota. The set-up in this study consisted in a main bioreactor containing the luminalassociated microbiota and a connected glass compartment with mucin beads to simulate the mucusassociated microbiota. The system was operated to simulate the colonic conditions of healthy human adults as described earlier (temperature 37°C, pH 6.3, retention time 24 h) 25,26 . The experiments were conducted in duplicate with fecal samples from two donors (one male and one female, ranging in age from 24 to 50 years, with no history of antibiotic or probiotic treatment 3 months prior the beginning of the study) 25 . Following fecal inoculation of the bioreactor, fermentations were conducted for a total duration of 9 days, including 1 day under fed batch and the following 8 days under semi-continuous mode. Samples were collected daily in the bioreactor.
In situ metagenomic measurement of growth rate by Growth Rate Index (GRiD) scores The GRiD software (v1.3) 27 was used to calculate the growth rate index from the metagenomic samples from Swedish wellness cohorts and fecal samples inoculated into bioreactor and fermented after 24 hours. Brie y, this software calculates the growth rate by mapping the metagenomics reads to microbial genomes and calculates the coverage ratio between the origin and terminus of replication as a proxy for the growth rate. Since GRiD is sensitive to the representativeness and quality of the genome used, we created a GRiD custom database representative to the gut microbiota, consisting only of high-quality draft genomes from the MGnify database 28 . First, we matched the msp gene clusters to the MGnify genomes using a BLASTN procedure, with a 95% identity threshold. Then we kept only the MGnify genomes passing these criteria: (1)  and KEGG orthology (KO) annotation of the gut catalogue. The KO pro le of each MSP were mapped into KBase metabolic model 30 as reference model to provide reaction pro les. Regarding the reaction pro les the context speci c GEMs were reconstructed and the functionality of the models was checked based provided biomass objective function and the gap lling was done using the COBRA toolbox and the reference model. To investigate the response of the IFS and OFS MSPs to environmental changes and calculate the perturbations, we used four different diets i.e. high protein-and bre-plant based diets and high-protein and bre omnivorous diets. The composition of the diet was converted to mmol/gDW*hour for the simulation in anaerobic situation and the growth rate for each model were predicted for each diet using constraint-based modelling. To check the dependence of the IFS and OFS species to the compounds as input or medium and autotrophy, we performed an essentiality analysis in which the inability of each MSP to synthesize the metabolites was simulated by closing the corresponding exchange reactions; decreased growth rate shows the dependence of the MSP to the metabolites for growth.
Gene richness and species associated with high and low gene richness Unsupervised clustering of co-conserved functions of gut microbiota We calculated Jaccard index among functional annotations to check how many species were sharing given a pair of functions together. We selected highly shared pairs of functions (Jaccard index >= 0.75) and merged into functional co-occurrence network using R igraph package 31 . Functional clusters within the network were identi ed by unsupervised community detection, short random work algorithm (cluster_walktrap function) 32,33 and identi ed singleton functions within the network as well. Among nonsingleton functional clusters, we selected representative functional clusters if functions of given functional clusters were found more than three species, thereby excluding functional clusters sparsely annotated over MSPs. Associated MSPs to functional clusters were chosen if given MSP covered more than 75% functions of given functional cluster.

Declarations
Data availability All primary HGMA data together with the newly sequenced data are freely available without restrictions in the public open access database (https://www.microbiomeatlas.org) that is part of the Human Protein Atlas program (https://www.proteinatlas.org).

Code availability
The R package used to perform modelling temporal changes of microbiome for in ow and out ow analysis together with functional clusters including unsupervised clustering of co-conserved functions of gut microbiota can be found at our GitHub repository link: https://github.com/sysbiomelab/mPackage.
The modeling of temporal changes can be applied directly to any sets of longitudinal microbiome data. The functional cluster analysis can be applied on gene counts and species abundances. The other pipeline scripts for analysis are also publicly shareable and available upon reasonable request from the corresponding author.  Square concerning the legal status of any country, territory, city or area or of its authorities, or concerning the delimitation of its frontiers or boundaries. This map has been provided by the authors. butyriciproducens and B. obeum), out ow species (OFS) (e.g. V. parvula and H. e uvia) or species that often change state between visits (E. coli, S. salivarius). b, we determined abundance changes between visits (t and t+1) and compared the total changes of IFS and OFS (ΔIFS and ΔOFS, respectively). In short, ΔIFS and ΔOFS were determined by differences of scaled total abundance between visits (i.e. Δ = Zt+1 -Zt) (See Methods). We observed negative correlation of abundance changes by visit between IFS species (ΔIFS) and OFS species populations (ΔOFS) (Spearman's ρ = -0.334, p-value = 4.6×10-8). c, we compared the mean abundance of OFS species in visits (μOFS) and abundance changes between visits (∆OFS).
The mean abundance of OFS species in visits was determined by the mean of scaled total abundance of Page 25/27 OFS species in two sequential visits (i.e. μOFS = ½ × (Zt + Zt+1)). The higher the mean OFS species abundance (μOFS ≥ 2), the more its abundance changes between visits either increasing (∆OFS > 2) or decreasing (∆OFS < -2). d, Individuals enriched with OFS species at certain visits (μOFS ≥ 2) had signi cantly lower intra-individual similarity (Wilcoxon one-sided test P = 0.0034). e, we observed that individuals with an increase in IFS species abundance were signi cantly positively correlated with an increase in richness (Spearman's ρ = 0.206, p-value = 9.0×10-4). We estimated physiological properties of IFS and OFS species by growth rate estimations by GRiD scores (f-i) and genome-scale metabolic modelling (j-k). We estimated GRiD scores of IFS and OFS species from (f) individuals highly enriched with IFS species and (g) individuals highly enriched with OFS species and observed higher GRiD scores of IFS species. In additional experiments of bioreactor fermentation of human faeces, we observed higher GRiD scores of IFS species after 24 hours, compared to original feces samples, whereas OFS species remained to be lower in GRiD scores at 24 hours. We evaluated our ndings in ow/out ow scores estimated from Swedish cohort of four visits with those from (l-m) additional two more visits, (n-o) American HPFS cohort and (p-q) Italian DINAMIC cohort. We compared in ow scores (l, n, p) and out ow scores (m, o, q) between different datasets. and negative associations (-) are marked on a heatmap (size proportional to signi cance). Increase of IFS species abundance was associated with increase of muscle mass, whereas increase of OFS species was associated with increase of BNP, a heart failure marker.

Figure 4
In ow and out ow were associated with disease signatures from pan-disease analysis, together with biodegradation of xenobiotics and host nutrient hijacking. a-b, we observed signi cant differences in in ow (a) and out ow (b) probability between common enriched and depleted MSPs in diseases (from Fig. 1h) (Wilcoxon one-sided tests p-values < 0.05). we observed that common depleted species were higher in in ow, whereas common enriched species were higher in out ow. c-d, radar plots showing the fraction of functional classes enriched in either (c) in ow or (d) out ow, tested by linear mixed effect models (adjusted p-value < 10-3). The in ows were enriched in core metabolism and out ow in accessory metabolism (e.g., BTEX contaminants). e, we identi ed functional clusters co-conserved across different species and distinguished functional clusters found in many species or few species (x-axis). Y-axis indicates the size of functional clusters. Here we identi ed two largest clusters enriched in pathobionts (CL-10, patho-cluster) and commensal bacteria (CL-12, comm-cluster), respectively. f, frequency of functional clusters associated with enriched/depleted MSPs in diseases: total frequency of enriched cohorts and depleted cohorts per functional cluster (|number of enriched cohorts| + |number of depleted cohorts|, y-axis) and subtracted frequency between enriched cohorts and depleted cohorts per functional cluster (|number of enriched cohorts| -|number of depleted cohorts|, x-axis). Functional clusters signi cantly associated with out ow species (hypergeometric tests P < 0.01) were mostly associated with depleted MSPs in diseases (points coloured purple). Points were plotted in a jittered setting to avoid overlaid points.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.