Metagenomic analysis of microbe-mediated vitamin metabolism in the human gut microbiome

Metagenomic abundance of vitamin biosynthetic pathways and transporters in healthy individuals from different countries

To evaluate if the abundances of microbial-mediated vitamin metabolic genes are dependent on the country of origin or not, linear and logistic regression analyses, without controlling for other confounding factors, was carried out. The result suggested an association between vitamin metabolic gene abundances with healthy samples originating from China (F-statistic, p < 0.01, as a country variable). However, the gene abundances from healthy samples of other countries, i.e., Denmark, Spain, and USA showed no significant association to the corresponding country variable. Extraction of DNA was performed in the Beijing Genome Institute (BGI) for the samples from China, Denmark, and Spain. Although none of these studies provided diet information for the samples, it would be reasonable to speculate that the observed differences could be a diet-dependent effect (based on Western and non-Western diets).

Henceforth, we compared the normalized abundances of eight B-vitamin and menaquinone related biosynthetic and transporter genes in the healthy subjects from these different populations. On comparison of healthy samples among four different countries, we found that individuals from China had significantly higher gene abundance (Mann–Whitney test, p < 0.01) for each vitamin type, than that from USA, Spain, and Denmark (Fig. 1 and Additional file 1: Table S1). For each vitamin, the normalized abundance of the biosynthetic genes ranged around similar range as shown in Fig. 1 for all the healthy individuals irrespective of their origin of country. Additionally, genes involved in the biosynthetic pathways were found to be significantly higher than their corresponding transport related pathways (Mann–Whitney test, p < 0.01). However, it is of caution that as the bacterial vitamin transporters could belong to a diverse class of protein families [15, 16] due to their evolutionary variability, functional annotation of these transporters through bioinformatic approach remains challenging [17, 18].

Metagenomic abundance of vitamin biosynthetic pathways and transporters in different states of health

Many functional features of gut microbiota have been shown to correlate with health and disease, and we therefore evaluated if vitamin metabolism by gut commensals is associated with the host health status. To do so, we performed similar analysis as described above between country-matched healthy controls and diseased subjects. The comparison for gene abundances included the following groups: healthy controls and the newly-onset IBD (inflammatory bowel disease) subjects from USA [14]; healthy controls and T2D (type 2 diabetes) subjects from China [13]; and healthy controls and IBD subjects under clinical remission from Spain [19]. The results suggested that gene abundances from T2D subjects in China (as a diagnosed variable) showed an association for vitamin metabolism (F-statistic, p < 0.01), while the gene abundances from IBD subjects from USA and Spain showed no significant association with disease status.

Comparing the gene abundances between the country-matched groups, we found that the T2D subjects from the Chinese cohort showed substantial and significant difference for all the vitamin types associated with biosynthetic enzymes and transporters (except for biotin and riboflavin) in comparison with the country-matched healthy controls (Mann–Whitney test, p < 0.01) (Fig. 2). On the other hand, IBD subjects (both onset and clinical remission ones) were not found to be significantly different concerning vitamin metabolism when compared to their respective country-matched healthy controls. To validate our analytical method, we tested and found that our findings were consistent with what has been reported earlier [20]. Henceforward, data for healthy subjects were used to get insight into how these genes map to bacteria from different phyla.

Phylogenetic origin of vitamin pathway genes

To compare the taxonomic profile of vitamin metabolic genes (biosynthesis and transporter genes) at phylum or species level, each of the mapped gene identifiers were used to retrieve the lineage information from UniProt database. In Fig. 3, we show how vitamin metabolic genes are prevalent across different phyla, which might suggest vitamin metabolism being as a part of core function in gut microbial metabolic modules. Additional file 2: Figure S1 shows a network of associations between species with mapped vitamin-related gene and its corresponding phylum. This figure reveals that the vitamin related pathways are mostly found in the phylum Firmicutes with about 49%, followed by Proteobacteria with about 19%, Bacteroidetes with about 14% and Actinobacteria with about 13%. Additional file 3: Figure S2 shows a network of associations between vitamin biosynthetic capabilities and phylum with species information. For every vitamin type, the number of species having prevalence of a biosynthetic pathway seemed to be higher than that of species with transporter mechanisms. Among the various vitamin transporters, species encoding thiamine transporter dominated in contrast to folate transporter being the lowest in number. However, from Additional file 3: Figure S2, it also becomes apparent that there are three types of species that differ with respect to functional characteristics related to vitamin metabolism, i.e. vitamin genes related to biosynthesis or transporters or both. Based on this, we categorized 171 microbial species as Producer (P), Consumer (C) and Producer-Consumer (Dual) for three specific vitamins (biotin, cobalamin, and thiamine) in the following analysis.

Network association between vitamin prototrophs and vitamin consumers

Given a hypothetical situation where there is a direct dependence between public-good-producing bacteria and its consumer, it is expected to observe a stronger association of their coexistence at the level of species abundance. Based on the presence or absence of genes encoding proteins required in vitamin biosynthetic pathways or transporter mechanisms or both, each microbial species was assigned as either vitamin prototrophs (producers) or vitamin consumers (consumers) or producer-consumer (dual) respectively for every vitamin type. Therefore, to estimate interactions between vitamin prototrophs and consumers, we performed a correlational analysis between species with all possible phenotypes. Of total 7021 possible pairs of species combinations (171 species (675 strains) considered), 2030 (i.e. 28.9%) and 590 (i.e. 8%) were significantly positive and negative correlated, respectively (Mann–Whitney test, p < 0.01). Among the positive correlations, with a correlation coefficient greater than 0.7, we observed that a pair of species originating from the same phylum could coexist, probably due to their evolutionary similarities. However, we also observed a similar level of coexistence for a pair of species originating from different phyla. For example, species from two different phyla, Lactobacillus acidophilus (Firmicutes) and Citrobacter koseri (Proteobacteria) had a high and positive correlation coefficient r = 0.81.

Based on the high prevalence of species with vitamin metabolic pathways among different phyla as derived from Fig. 3, we focused our next correlational analysis (using similar approached as described above) to the most prevalent three B-vitamins among microbial species, i.e. biotin, cobalamin, and thiamine involving both biosynthetic and transporter related genes respectively to estimate interactions between these specific vitamin prototrophs and consumers. Contrary to our hypothesis, we found that species with dual functional traits (i.e. containing genes for both biosynthesis and transport of the same vitamin) contributed to a notable fraction of positive correlations among the possible phenotypes (Additional file 4: Figure S3). On the other hand, the co-occurrence of species with single phenotype (i.e., either producer/P or consumer/C) and dual phenotypes belonged to both positive and negative correlations. These observations hint that the species with dual phenotypes could drive stronger relationships than the ones with either limited or less functional traits. It is well-known, that communities of species with multiple functional features respond well to changes in the environment, thereby influencing the net outcome of a community assembly [21–23]

Genome annotation coverage of biosynthetic pathways

Next, we sought to visualize the presence and absence of vitamin pathway related reactions in each cohort to gain insights into how conserved each pathway is, across different populations. As depicted in Additional file 5: Figure S4, riboflavin seems to be well represented and conserved among all species across all the four cohorts regardless of the health status. This conservation could also be due to the limited number of reactions involved in riboflavin pathway compared to other pathways. Contrarily, menaquinone or vitamin K₂, has the least number of species with its biosynthesis pathway present. Perhaps, human microbiome might have evolved to have a reduced dependence of menaquinone production from the gut microbiota as proposed earlier [24]. In case of folate, the majority of the reactions involved in its biosynthesis were observed to be present except for four reactions which concerned interconversions of tetrahydrofolate, dihydrofolate, and folate using NADH or NADPH as cofactors and one ATP-forming reaction carried out by folylpolyglutamate synthase enzyme. Biotin and thiamine appeared to be vitamins with a major part of their biosynthetic reactions with their annotations missing. These are related to fatty acid metabolism and the upstream part of the biotin synthesis pathway provides a long-chain acyl as the starting point for the biosynthesis. In case of thiamine, three phosphatase reactions involving interconversion of thiamine monophosphate, −diphosphate, −triphosphate and thiamine were absent. For the rest of the vitamins, the incompleteness of their biosynthesis pathway was observed for many species, which is in line with previous studies as it has been suggested that vitamin biosynthesis could be carried out as complementary tasks between several bacteria, each harboring a part of the pathway [25].

Furthermore, for the American and Chinese cohort, we investigated the presence and absence of reactions in the biosynthesis pathway of each vitamin considering the abundance of bacterial species. As depicted in Fig. 4a, it is noticeable that species with biosynthetic reactions of riboflavin, pantothenate, and pyridoxine have relatively higher coverage when compared to the other vitamin types. It is also to be noted that species from genera Bacteroides, namely B. uniformis, B. stercoris, B. vulgatus, B. fragilis, and B. caccae are among the most abundant species with relatively high coverage of biosynthetic capabilities for all vitamins. Although the overall coverage of vitamin biosynthetic pathways is rather high, the same trend is not observable in case of biotin, cobalamin, and niacin which have incomplete pathways in many cases.

For the Chinese cohort (Fig. 4b), we noticed that the abundant species profile differs from that of the American cohort, possibly reflecting a population-specific (which is mostly due to country variation) gut microbiota. It is again evident from this that Bacteroides species have most of the vitamin biosynthesis pathways covered except for cobalamin due to differences at the species level. Besides, in the Chinese population, there are many species from Firmicutes amongst the abundant ones, which suggest that based on the Bacteroidetes/Firmicutes ratio there may be a change in abundance of vitamin related genes.

Expression level of vitamin biosynthetic enzymes and vitamin-dependent enzymes in an IBD cohort

Based on the availability of metatranscriptomics dataset from the American cohort analyzed above, we sought to investigate the mRNA sequences of microbe-mediated vitamin related enzyme as computed by HUMAnN2 [26]. Functional profiling of mRNA normalized features (i.e., vitamin biosynthesis) in this cohort reveals that there is a low abundance of cobalamin and thiamine biosynthesis related enzymatic genes expressed in the IBD samples compared to healthy controls (Mann–Whitney test, p < 0.01) (Fig. 5a). These findings were consistent with a similar pattern observed in a previous study where B-vitamins were measured to be lower in plasma samples of IBD subjects [27–29].

However, the expression level of biotin biosynthetic enzymes was found to be significantly higher in the IBD subjects compared to healthy controls (Mann–Whitney test, p < 0.01). To examine this, we inspected whether the abundance of biotin synthesizers were significantly abundant or not in the IBD subjects than that of the healthy controls (especially, for the coordinated set of samples where the metagenomics and metatranscriptomics were carried out). Interestingly, we found that species with only biosynthetic pathways for biotin were significantly more abundant in the IBD subjects whereas species with biotin transporters were significantly less abundant in the IBD subjects compared to healthy controls (Additional file 6: Figure S5).

These observations once again highlight the significance of the microbial contribution to human health status. Furthermore, we analyzed the expression level of these vitamin-dependent enzymes to discern the connection between their synthesis and usage. We found that the abundance of biotin-dependent enzymes had an inreverse pattern to that of their biosynthesis level, whereas, thiamine utilization potential as well as its biosynthetic potential were less abundant in IBD subjects when compared to the healthy controls (Fig. 5b) (Mann–Whitney test, p < 0.01). Variation in vitamin synthesis and utilization level among these three vitamins is likely to reflect the complexity of the gut microbiota towards resource allocation and reuse of enzymes.

Cohort description

Metagenomic sequences subjected to Illumina paired end sequencing were obtained from European Nucleotide Archive at EMBL-EBI (https://www.ebi.ac.uk/ena) under the accession number PRJEB2054, PRJNA422434, PRJNA275349, and PRJNA389280.

From PRJEB2054, 14 and 25 samples were of healthy controls and IBD subjects from Spain respectively [19]. Healthy controls were recruited among family relatives of IBD patients; antibiotic treatment for at least 4 weeks before fecal sample collection was excluded. IBD subjects were in clinical remission for at least 3 months, and had stable maintenance therapy with mesalazine or azathioprine [33]. 15 and 70 samples originating from Denmark were of low-gene count and high-gene count respectively [34].

From PRJNA422434, 187 healthy controls and 172 subjects who were diagnosed with Type 2 Diabetes were analyzed. These fecal samples were sampled from Shenzhen Second People’s Hospital, Peking University Shenzhen Hospital and Medical Research Center of Guangdong General Hospital, individuals living in the south of China. Detailed clinical information for the samples is available in the supplemental information of the original article [13].

From PRJNA389280, 55 healthy controls and 223 IBD subjects were analyzed. None of the IBD subjects had history of terminal ileal resection. The samples were collected and recruited from Massachusetts General Hospital, Cedars-Sinai Medical Center, and Emory University and Cincinnati Children’s Hospital [35]. Detailed information on the criteria considered for the sample inclusion or exclusion on IBD subjects is available in the original article [14].

Functional mapping and analysis

Vitamin related pathway abundances were calculated using an open-sourced, stand-alone functional analysis pipeline FMAP [36]. Metagenomic read sequences were mapped to the reference database using DIAMOND [37] as mapping mode. Unique hits that passed through the default parameters (i.e., e-value >1e-10 and sequence identity > 80%) of the pipeline FMAP were considered for further analysis. From the BLAST output files, the subject sequence identifiers were used for identification of taxonomic lineage and KO analysis. Gene sequences as determined by UniRef were mapped to Kyoto Encyclopedia of Genes and Genomes (KEGG) Orthologies (KOs) database [38, 39] and were binned into functional categories as listed in Additional file 7: Table S1-S9. Taxonomic lineage of the mapped hits were retrieved from the UniProt website (dated on October 2018). For every vitamin type, total normalized abundance, was computed as number of reads that have been aligned to a gene divided by total number of read counts in a sample followed by summation of all the genes constituting the pathway. Or in other words, the summation of the normalized abundances of each vitamin’s KO in each sample was defined as the total normalized abundance for each vitamin pathway.

Relative abundance of the species was calculated using MEDUSA [40] for the MetaHIT and T2D cohort. Data table with relative abundance of species and metatranscriptomics data for the HMP2 cohort were obtained from their publicly available project website [14]. Using the species abundance values, spearman correlation was applied to estimate correlations in determining the relationship between different microbial species. Pairs of species with correlation coefficient of greater than or equal to 0.4 were considered as positive correlations. Pairs of species with correlation coefficient of less than zero were considered as negative correlation. Based on the taxonomic information extracted for each of the mapped protein identifiers from the UniProt, microbial species encoding proteins required for biosynthesis of a certain vitamin but lacking transporters for it were assigned as vitamin prototrophs (producers or ‘P’). The species having the transporters but lacking the biosynthetic genes for a certain vitamin were assumed as vitamin auxotrophs (consumers or ‘C’). As well, the species containing genes for both the biosynthesis and the transport of the vitamin were assigned as dual (producer-consumer or ‘PC’).

For identification of vitamin-dependent enzymes, their enzyme commission numbers were obtained from ExPASy (listed in Additional file 7: Table S10).

Metabolic reconstructions

In order to investigate the coverage of vitamin biosynthesis pathways, previously published repository of human gut bacterial genome-scale metabolic models were deployed [41]. These models have a GPR (gene-protein-reaction) association component which were used to translate gene profiles to reactions through logical associations. From this repository, vitamin biosynthesis profile of the overlapping species was looked in each cohort. To do so, all the KO groups related to vitamin biosynthetic reactions as described in the previous section were extracted. KEGG reactions were mapped to Kbase [42] nomenclature systems to be consistent with the reaction identifiers in the metabolic models. Following identification, a list of biochemical interconversions for each vitamin category, namely, biotin, cobalamin, folate, niacin, pantothenate, pyridoxine, riboflavin, thiamine and menaquinone were complied. Subsequently, the presence and absence of these reactions in each cohort were analyzed to see how conserved each pathway is across different populations. Furthermore, for the American and Chinese cohort, we investigated the presence and absence of the reactions in biosynthesis pathway of each vitamin considering the abundance of bacterial species. In each community, a union list of the present species was constructed and sorted according to the relative abundance of the species in the healthy cohort. Species unique to the diseased population were then added to the end of the list.

Statistical analyses

All the statistical analysis and graphics were performed using R version 3.5.1 [43]. Country matched samples were compared into two groups for differential KO abundance testing between healthy controls and diseased subjects to avoid cross country-effects. Statistical analyses that involved two groups and more than two groups was performed by Mann-Whitney-Wilcoxon tests and Kruskal Wallis tests respectively. Bonferroni’s correction was applied to the significance level. For multivariable linear regression, categorical metadata and vitamin type abundances, samples were pooled into classifications (male/female, healthy/diseased, specific nationality/rest, etc) and significance was identified using Fisher’s exact test with multiple testing correction of P values.