Dataset collection and meta-analysis workflow
In this meta-analysis, we collected 1711 samples from 16 metagenomics studies spanning 12 host phenotypes, including 919 nonhealthy and 792 healthy human individuals (Figure 1A, Additional file 2: Table S1). Notably, subjects were categorized into “healthy” or “nonhealthy” according to the original studies, while overweight and obesity alone individual were not classified as nonhealthy group. Expectedly, we observed beta-diversity differences based on species-level fecal microbial abundance between healthy and nonhealthy groups using Bray-Curtis dissimilarity (R=0.002, p<0.001, Adonis, Figure 1B, 95% confidence regions), consistent with the past results in most studies.
Subsequently, we sought to explore the disease-associated genetic variation in gut microbiota by mapping metagenomics reads against a comprehensive set of gut microbial reference genomes. If we identified any SNVs from a microbial reference genome, this “mutated genome” or SNV-carrying microbial strain will be further analyzed to characterize its SNV profile and test if it can associate with disease status (Additional file 2: Table S2). We identified a total of non-redundant 2740 mutated genomes in 1711 samples. Here, we calculated the distribution of average coverage and breadth of these mutated genomes with prevalence greater than 10% (233, 8.50%), as meaningful comparison of genomic variation between samples requires both breadth and depth of sequencing for each genome [6]. In this study, the sequencing coverage of each SNV was required to be 5, minimally, which could cover most genomes (average coverage of 158 genomes more than 5, 67.81%) (Additional file 1: Figure S1). Interestingly, we next observed the beta-diversity difference based on the composition of SNV-carrying strains (N=2740) between healthy and nonhealthy groups (R=0.003, p<0.001, Adonis, Figure 1C, 95% confidence region) based on Euclidean dissimilarity matrix. This suggested that host medical conditions can associate with not only population-level abundance changes but also genetic compositional changes in the gut microbiota. Next, we focused on specific strains, genes and SNVs, showing the association between gut microbial SNVs and host health status.
The healthy gut microbiome harbors a wider range of SNV-carrying strains
According to InStrain pipeline for analysis of co-occurring genome populations (“Methods”), we obtained the comprehensive SNV profiles of gut microbiota for both healthy and nonhealthy cohorts. Of note, in this study, the number of strains and SNV count was normalized based on metagenomic sequencing depth (“Methods”, Additional file 1: Figure S2).
We first explored whether any difference in the richness of mutated strains, and the total number of SNVs between healthy and diseased states. Firstly, we observed that healthy subjects (19.48±7.38) had more mutated strains than diseased subjects (17.39±6.84) (Figure 2A, Additional file 2: Table S3, p<0.001), suggesting a higher strain-level diversity in the healthy gut microbiota. Respectively, such a difference pattern was also found in a total of six independent cohorts corresponding to six diseases, including AS (16.92±5.64 vs. 21.93±7.23, p<0.001, average±std), CD (17.33±6.87 vs. 24.76±5.77, p<0.001), GD (14.78±3.75 vs. 17.61±3.37, p<0.001) and LC (19.51±6.01 vs. 23.77±6.51, p<0.001) and SF (17.56±3.81 vs. 25.51±3.0, p=0.016). In contrast, CRC individuals had a higher strain-level diversity than healthy control (22.26±6.80 vs. 20.73±7.11, p=0.025). There were no significant differences in other cohorts, including BC (11.10±2.66 vs. 11.83±2.83), Cov19 (14.29±7.41 vs. 16.28±10.84), PCOS (15.32±3.87 vs.15.6±14.67), SCZ (11.56±3.71 vs.11.29±3.45), T2D (19.13±5.61 vs. 19.78±5.18) and UC (22.91±8.62 vs. 24.76±5.77).
Most studies have reported differences in the species-level alpha diversity of healthy and nonhealthy cohorts. Here, our integrated analysis showed that alpha diversity (using either Shannon or Simpson index) is different between host disease states (Figure 2B, p<0.001, Wilcoxon rank-sum test). Intriguingly, the richness of mutated genomes strongly correlates with Shannon index (R=0.652, p<0.001), or Simpson (R=0.596, p<0.001) (Figure 2C. One possibility is that lower intestinal microbial diversity leads to a smaller number of mutated genomes. To address these concerns, we matched healthy and nonhealthy patients with alpha diversity within 0.01 (Shannon index) or 0.001 (Simpson index) differences to compare their normalized number of mutated genomes (92.87% and 92.69% samples were included), and Shannon and Simpson index were 0.01 or 0.001 higher in nonhealthy group than in healthy. In such manner, we found that when the alpha diversity is almost the same, the number of mutated genomes is different (Wilcoxon rank-signed test, p<0.001 for samples with equalized Shannon index, and p=0.023 for samples with equalized Simpson index). Therefore, in most cases, we believe that healthy individuals have higher strain-level diversity than nonhealthy individuals (Figure 2D, Additional file 2: Table S4). At the same time, we also matched individuals with < 0.01 difference in the normalized number of mutated genomes, and compare their alpha diversity (27.47% samples were included) between host groups. Interestingly, we still observed the disease-related difference in Simpson index at the species level (p=0.032), where healthy group had less microbial species (Figure 2D, Additional file 2: Table S4). Therefore, healthy subjects can reach equal strain-level diversity with nonhealthy subjects when they had less microbial species in the gut microbial community. Overall, healthy individuals tend to have higher strain-level diversity than most nonhealthy individuals, expect CRC patients.
Next, we estimated the SNV number per Mkb sequencing data for all mutated genomes in both healthy and nonhealthy gut microbiota, and demonstrated healthy individuals have more SNVs in total than nonhealthy individuals (23.44±9.93 vs. 25.20±10.11, Figure 2E, p <0.001). Respectively, we validated such a difference in a total of five diseases/cohorts, including AS (212.85±67.13 vs. 264.80±65.15, p<0.001), CD (211.48±81.11 vs. 301.02±61.23, p<0.001), GD (213.46±54.13 vs. 241.71±56.14, p=0.012) and UC (257.86±107.20 vs. 301.02±61.23, p=0.002). By contrast, CRC individuals have more SNVs than healthy control (260.39±75.34 vs. 241.52±71.09, p=0.045). No significant differences were found in other cohorts, including BC (148.22±40.47 vs. 156.71±49.76), Cov19 (190.14±125.20 vs. 260.43±238.00), LC (270.74±77.99 vs. 289.17±64.32), PCOS (214.99±73.86 vs. 216.58±65.64), SCZ (148.39±62.71 vs. 135.49±58.73), SF (604.81±280.69 vs. 781.80±83.49) and T2D (371.49±78.86 vs. 339.52±79.23). Given the genetic compositional difference found at such a higher level between disease states, SNV profiles on specific strains/mutated genomes need to be revealed further.
Strain-level diversity associated with multiple host health status
To further explore if disease states can also associate with the SNV profiles for individual microbial strains and other key evolutionary patterns for gut microbiota, we comprehensively analyzed all 75 strains with a mutation frequency of more than 30% in the 1711-member population (Figure 3), and the corresponding genome ids of these strains can be found in Additional file 2: Table S2. Firstly, we compared the prevalence of these strains (Additional file 2: Table S2), and the top-four most prevalent strains in both healthy and unhealthy groups were: Bacteroides dorei (non-healthy: 87.16%, healthy: 91.92%, all: 89.36%), Bacteroides uniformis (non-healthy: 76.50%, healthy: 86.87%, all: 81.30%), Faecalibacterium prausnitzii G (non-healthy: 65.18%, healthy: 77.90%, all: 71.07%) and Parabacteroides distasonis (non-healthy: 69.42%, healthy: 71.21%, all: 70.25%). 46 strains were more prevalent (46/75, 61.33%) in the healthy group, while 5 strains (5/75, 6.67%) were more prevalent in the nonhealthy group. At the species level, many strains from Bacteroides consistently have a high frequency of SNVs occurred in the healthy gut. At the strain level, Faecalibacterium prausnitzii G, Faecalibacterium prausnitzii D, Faecalibacterium prausnitzii C, Faecalibacterium prausnitzii K and Faecalibacterium prausnitzii E also presented diverse patterns between healthy and nonhealthy population. Overall, the above results also confirm that the healthy group will have a wide range of mutations in the genome (Figure 2A).
We next investigated if SNV rate (“Methods”, Eq 1) and microdiversity of these gut microbial strains can associate with disease states. “SNV rate” is defined as, the relative frequency of SNVs on a given genome that can compare with other genomes. Microdiversity is a measurement of intra-population genetic diversity for a microbial strain in the gut microbiota (Additional file 1: Figure S3). Clearly, the SNV rate and microdiveristy substantially varied among microbial species even strains (e.g., Faecalibacterium prausnitzii). We identified 15 strains (out of 75 strains, 20%) that had a SNV-rate difference between disease states. Among these 15 strains, 8 strains (8/75, 10.67%) were higher in the nonhealthy group and 7 (7/75, 9.33%) strains were higher in the healthy group. Remarkably, eight out of 15 strains (8/15, 53.33%) were mainly from Lachnospiraceae, including Blautia wexlerae, Agathobacter rectalis, Anaerostipes hadrus, Lachnospira eligens, Blautia sp900066165, Faecalicatena torques, Roseburia intestinalis, KLE1615 sp900066985, a typical group of bacteria that produce SCFA (i.e., acetate and butyrate). In contrast, Bacteroides were the main genus (6 strains, 6/20, 30%) with microdiversity difference between disease states, including Bacteroides dorei, Bacteroides uniformis, Bacteroides xylanisolvens, Bacteroides sartorii, Bacteroides thetaiotaomicron and Bacteroides massiliensis, and only three strains (3/75, 4%, 2 Bifidobacterium) were enriched in the nonhealthy group while 17 (17/75, 22.67%) were enriched in the healthy group. Overall, these results confirmed that healthy individuals have more SNV in their gut microbial genomes (Figure 2E).
Universal base mutation bias associated with host health states
Next, we narrowed down to the specific pattern of mutation types at the single SNV level (i.e., nucleotide diversity) since the SNV rate and microdiversity of 75 gut strains were associated with disease states. Firstly, we profiled six mutation types of all 1771 samples (Additional file 2: Table S5), and compared alpha diversity and beta diversity based on mutation-type profiles across samples. Intriguingly, Shannon diversity for the mutation types (p=0.03) had significant differences between two groups (Figure 4A). In addition, the mutation-type composition of 75 gut microbial strains was primarily clustered by host disease states (R2=0.004, p<0.001, PCoA1=50.72%, PCoA2=12.36%) (Figure 4B). Specifically, we found that two mutation types (A>G|T>C, p=0.006 and A|T>T>A, p=0.011, Wilcoxon rank-sum test) were enriched in the healthy group while one in the nonhealthy group (C>G|G>C, p<0.001, Wilcoxon rank-sum test) (Figure 4C, Additional file 2: Table S5). We also showed and compared the mutation types between host states for each strain (Additional file 2: Table S6). Overall, we observed the diverse distribution patterns of mutation types among microbial strains (Figure 4D). C>G|G>C in many strains tend to be significantly higher in nonhealthy people (31/75, 41.33%), while such mutation type in only three strains were higher in the healthy group (3/75, 4%). In contrast, the A>G|T>C mutant type was more common in multiple strains in the healthy group (16/75, 21.33%).
We next compared the overall percentage of adaptive mutations occurring in the gene coding region for all 75 strains, and no differences between disease states were found (89.84% vs. 89.75%, p=0.1, Figure 4E). However, at the strain level, we observed mutations in coding regions of 32 strains were significantly higher in the nonhealthy group than in the healthy group (32/75, 42.67%), while only three strains (Bacteroides sp.) were higher in the healthy group (3/75, 4%). This suggested that disease states can change SNV patterns in coding regions of gut microbial strains.
Strain-level codon mutation bias in SCFA-production-involved genes
We next focused on key functional genes of specific gut strains and assessed potential functional changes in gut microbiota induced by microbial adaptive evolution under disease conditions. SCFA is widely considered to be closely related to human health, and SCFAs abundance deficiency has been observed in a variety of diseases, and the genetic evolution process of SCFAs related genes has not been reported so far. We then first counted the enzymes encoded by these strains that are related to SCFAs production, including acetate kinase (ack, E.C. 2.7.2.1), Propionyl-CoA:succinate CoA transferase (scpC, E.C. 2.8.3.-) and Butyrate kinase (buk, E.C. 2.7.7.7), acetate CoA-transferase YdiF (ydiF, E.C. 2.8.3.8), acetate CoA-transferase subunit alpha/beta (atoD/A, E.C. 2.8.3.8), butyrate--acetoacetate CoA-transferase subunit A/B (ctfA/B, E.C.2.8.3.9) and butanoate coenzyme A-transferase (E.C. 2.8.3.-, BcAt) (Figure 5A).
We next presented the SNV profile of 75 strains encoding the above SCFAs-related enzymes (Figure 5B). We observed that almost all strains encode ack gene (74/75, 98.67%), involved in the acetate production. 74.67% strains have at least one metabolic pathway associated with butyrate production, and six strains have at least two pathways. In addition, Bacteroides spp. are the main producers of propionate (16/21, 76.19%). Then, we compared whether the pN/pS values of these genes were different between healthy and nonhealthy cohorts to assess the size of intestinal selection pressure. Specifically, the pN/pS value of ack gene of 11 strains was different between host groups, suggesting this gene can be under positive selection due to medical conditions. Its pN/pS ratio in four strains was higher in the healthy group, while seven strains were higher in the nonhealthy group. Differential intestinal selection pressures related to diseases drove the production of butyrate by related genes in 10 strains. Only scpC gene of Bacteroides stercoris had a difference in pN/pS, and the nonhealthy group was higher than the healthy group. Overall, the adaptive evolutionary patterns of SCFAs-associated genes in different strains are diverse due to host health status.
Next, we explored if and how these SNVs specifically modulated the production of SCFAs. Here, we compared the frequency that a SNV on SCFAs-related gene that causes codon to become termination codon or the start codon was inactivated (Figure 5C). Intriguingly, those SNVs were mainly found in four strains (Faecalibacterium prausnitzii C, Faecalibacterium prausnitzii D, Bifidobacterium pseudocatenulatum and Bacteroides stercoris). The frequency of SNV-derived codon changes in the nonhealthy group was significantly higher than that in the healthy group due to the disease-related selection pressure (Additional file 2: Table S7). Respectively, two genes of Faecalibacterium prausnitzii C possessed codon mutation bias affecting initiation and termination codon, including ack (8.98% vs. 4.37%, p=0.007) and BcAt (2.20% vs. 0.69%, p=0.04). In addition, codon mutation bias also appeared in ack gene of Faecalibacterium prausnitzii D (1.08% vs. 0%, p=0.023), ack gene of Bifidobacterium pseudocatenulatum (1.68% vs. 0%, p=0.039) and scpC gene of Bacteroides stercoris (5.03% vs. 1.86%, p=0.007). Therefore, the adaptive mutations occurred in the above strains in the nonhealthy group would block in the expression of SCFAs-production genes. Although the SCFAs genes of some 6 strains were subjected to greater intestinal selection pressure in the healthy group, it did not cause an increase in the above codon mutation types (Additional file 1: Figure S4). We also compared the relative abundance of the four strains aforementioned between host states and found that the three strains were more abundant in the healthy group, including Faecalibacterium prausnitzii C (0.56% vs. 0.61%, p<0.001), Faecalibacterium prausnitzii D (0.51% vs. 0.43%, p<0.001) and Bacteroides stercoris (2.43% vs. 2.56%, p<0.007). Therefore, we concluded that not only did the nonhealthy individual have fewer SCFAs producers in the gut microbiota, but also these producers’ gene expression was also blocked due to adaptive mutations. Finally, we reported potentially common codon mutations in gut microbial strains (prevalent in > 3 host individuals, Figure 5D, Additional file 2: Table S7) affecting SCFA production, including three amino acids of Faecalibacterium prausnitzii C ack gene, one amino acid of Faecalibacterium prausnitzii C BcAt gene and two amino acids of Bacteroides stercoris scpC gene.
GMHI based on SNVs profiles indicates host health status
Gut Microbiome Health Index (GMHI) was proposed to develop for predicting disease presence (or absence) using species-level taxonomic abundance profiles derived from metagenomic sequencing data. Overall, microbial abundance essentially reflects the ecological changes related to disease states. Given our findings above, we hypothesized that SNV profile can also predict the host states, reflecting the degree and extent by which gut microbes adapt under multiple disease conditions.
Firstly, to test and validate how well the abundance-based GMHI can predict the host states, sequence abundance profile of the mutated genomes was used to construct GMHI (abundance-based GMHI, “Methods” and Additional file 3: Table S8). The relative sequence abundance of a total of non-redundant 2740 strains was calculated for 1711 metagenomics samples in the discovery cohorts and 446 samples in the validation cohort. Notably, we selected microbial strains differentially abundant in DNA sequence between diseased states as the biomarkers for calculating GMHI, which was slightly different from previous study [27] (“Methods” and “Code availability”). Among 471 selected strain-level markers, 222 were enriched in the nonhealthy group, while 249 were enriched in the healthy group. We found that the overall prediction accuracy of the abundance-based GMHI was 73.97%, p=9.98e-66, while this GMHI showed a significant difference between healthy and nonhealthy groups in nine cohorts (Figure 6A). The result was similar to the reported accuracy of the original GMHI method (~70%) based on the taxonomic abundance of gut microbes at the species level [27].
We next explored if SNV rate, representing the scale of genetic variation for each microbial strain (Eq 1), can be used as an independent indicator to predict host disease states. Therefore, we calculated the SNV rate of each of non-redundant 2740 strains in each sample, which formed a new SNV-rate feature table as compared to conventionally used abundance-based feature table. Next, we identified microbial strains that have a significant difference in the SNV rate between host states and used their SNV rate to construct a new GMHI for classifying disease states (“Methods”, Additional file 4: Table S9). Among 470 selected markers, 214 were enriched in the nonhealthy group, while 256 were enriched in the healthy group. With the receiver operating characteristic curve (ROC), we further showed high prediction accuracy of this SNV-rate GMHI (AUROC=74.23%, Wilcoxon rank-sum test, p=4.12e-67). Its performance was slightly better than the strain-level abundance-based GMHI. Furthermore, we found this GMHI can distinguish host states in 10 independent cohorts (Figure 6B).
We’ve shown that the mutation frequency/bias of SCFAs genes in gut microbes can potentially distinguish health states. Therefore, we next explored the feasibility of constructing GMHI based on the SNV profile in SCFAs genes of 75 microbial gut strains. We collected 185 genes from nine gene families (ack, buk, ydiF, atoD, atoA, ctfA, ctfB, BcAt, scpC) from 75 strains that carried adaptive SNVs. We first formed a count table for the SNVs (presence and absence) identified in the 185 genes. Then we aggregated the SNV count to the genome level, and calculate the relative frequency of SNV for SCFAs genes in each genome in each sample (“Methods”, Additional file 5: Table S10). Among 48 selected markers, 11 were enriched in the nonhealthy group, while 37 were enriched in the healthy group. With the selected microbial SNV biomarkers, this SCFA GMHI showed that fairly good prediction accuracy (AUROC=65.20%, Wilcoxon rank-sum test, p=1.07e-26). Among 12 cohorts in this meta-analysis, the SCFA-based GMHI can be validated in six cohorts (Figure 6C). Although this accuracy is not perfect, it strongly suggested that genetic variability in the SCFAs genes alone can explain, at least partially, the microbiome difference between host health conditions.
External validation of functional changes in SCFA production induced by SNVs and the predictive power of GMHIs
To validate the above results, we further included three independent cohorts with a total of 446 samples, including 244 non-healthy and 202 healthy individuals, related to Atherosclerotic cardiovascular disease (ACVD) and Tuberculosis (TB) (Additional file 6: Table S11). The differences in relative abundance of the four strains changed compared with Figure 5 due to the different cohorts (Figure 7A). We observed that the relative abundances of Faecalibacterium prausnitzii C (p=0.37) was no longer significantly higher in the healthy group than in the non-healthy group, while the abundances of Faecalibacterium prausnitzii D (p<0.001), Bacteroides stercoris (p=0.042) and Bifidobacterium pseudocatenulatum (p<0.001) were higher in the healthy group. Interestingly, we still observed differential patterns of codon mutation bias in ack gene of Faecalibacterium prausnitzii C (p=0.023) and scpC gene of Bacteroides stercoris (p=0.025), however, other codon biases have not been proven (Additional file 1: Figure S5). We speculated that the other three genes were not verified due to the low probability of codon mutation bias (initiator and terminator), which can also be confirmed from Figure 5. Overall, we demonstrated that SCFAs metabolic genes on specific gut microbes exhibit codon mutation bias (initiator and terminator) in hosts with different health states.
Next, the ACVD and TB cohorts were further validated for GMHI based on SNV profiles (Additional file 6: Table S12-14). Overall, the three different GMHIs showed acceptable prediction accuracy, respectively. The validation accuracy (i.e., AUROC) of abundance-based GMHI was 71.50% (Wilcoxon rank-sum test, p=5.2e-15), that of SNV-rate GMHI was 69.28% (Wilcoxon rank-sum test, p=2.3e-12) and that of SCFAs-SNV GMHI was 63.57% (Wilcoxon rank-sum test, p=8e-7). Furthermore, these GMHI can distinguish ACVD or TB from health. Collectively, we demonstrated that health status changes can be distinguished by genetic changes in gut microbiome, and highlighted the importance of SNV-induced functional changes related to SCFA production which underlie the pathogenesis of inflammatory bowel diseases and many diseases that were not interrogated before.