Shotgun Metagenomics of 250 Adult Twins Reveals Genetic and Environmental Impacts on the Gut Microbiome

Abstract

The gut microbiota has been typically viewed as an environmental factor for human health. Twins are well suited for investigating the concordance of their gut microbiomes and decomposing genetic and environmental influences. However, existing twin studies utilizing metagenomic shotgun sequencing have included only a few samples. Here, we sequenced fecal samples from 250 adult twins in the TwinsUK registry and constructed a comprehensive gut microbial reference gene catalog. We demonstrate heritability of many microbial taxa and functional modules in the gut microbiome, including those associated with diseases. Moreover, we identified 8 million SNPs in the gut microbiome and observe a high similarity in microbiome SNPs between twins that slowly decreases after decades of living apart. The results shed new light on the genetic and environmental influences on the composition and function of the gut microbiome that could relate to risk of complex diseases.

Keywords: SNPs; heritability; metabolic diseases; microbiome.

Figure 1. Representation of TwinsUK Samples by the Gene Catalogs

High-quality non-human metagenomic reads uniquely (green, red) or commonly (blue) aligned to genes from the 250 twins gene catalog and the IGC (Table S1B). The average alignment ratio to each part is shown in the middle. The updated TwinsUK reference gene catalog (1,517 samples + 511 genomes, Table S1C) allows on average 80.21% mapping of the reads (unique + common). This is close to saturation because the percentage of gene-coding regions in all prokaryotic genomes is ~87%, and an estimated 7.25% of sequencing reads with an average length of 77 bp could not be mapped reliably as they only partially overlapped with genes (Li et al., 2014). Age and location of the samples are shown for reference.

Figure 2. Effect of Current Location on the Gut Microbiome

(A) Distribution of the subjects among geographic locations (Table S1A). The longitude and latitude of each location were hierarchically clustered to yield four regions (geo-clusters), represented by different colors, except for one sample from Ireland (not shown in B and C). The size of each circle scales with the number of subjects (samples) from that county. (B) Venn diagram for the number of the twin pairs in the same geo-cluster and in two different geo-clusters. (C) Bray-Curtis distance of the gut microbial gene profile between any two samples in the same region (intra-clusters), or in different regions (inter-clusters). Plotted are interquartile ranges (IQRs; boxes), medians (dark lines in the boxes), the lowest and highest values within 1.5 times IQR from the first and third quartiles (whiskers above and below the boxes), and outliers beyond the whiskers (circles). p = 1.99e-18, Wilcoxon rank-sum test. (D) Bray-Curtis distance of the gut microbial gene profile between paired twins in the same or different geo-clusters. p = 0.0363, one-tailed Wilcoxon rank-sum test. As the volunteer in Ireland has a twin sister in geo-cluster 2 (Table S1A), this pair is included in this panel.

Figure 3. Heritability of Gut Microbial Taxa

(A) A phylogenetic tree was drawn for the 90 genera seen in at least 50% of the samples. The heritability (A component in the ACE model) was plotted as a bar for each genus. Outer circle, green, ICC MZ < DZ; pink, ICC MZ > DZ; red, ICC MZ > DZ and p < 0.1 between ACE and CE models. Middle circle, colored according to phyla; inner circle, light to dark blue according to mean relative abundance of each genus. Genera that contained less than ten genes in 97% of the samples were not plotted. More detailed data are available in Table S3B for genera and Table S3A for phyla. (B) Heritability of mOTUs. The rank of relative abundance difference between each MZ or DZ twin pair, normalized to be between zero and 1, was shown as boxplots for Bacteroides, Bifidobacterium, Dorea, butyrate-producing bacterium, and unnamed mOTUs. Class information for all the plotted mOTUs was shown as colored bar on the left. Heritability of these mOTUs according to the ACE model were plotted to the right and color coded: pink, ICC MZ > DZ; red, ICC MZ > DZ and p < 0.1 between ACE and CE models. Detailed results for all mOTUs with more than 50% occurrence are available in Table S3C.

Figure 4. Heritability of Select Functions

(A–D) Heritability of KOs in the biosynthesis of branched chain amino acids (A), arginine and proline metabolism (B), phenylalanine, tyrosine, tryptophan, and histidine biosynthesis (C), and riboflavin metabolism (D). Pink, filtered KOs and ICC MZ > ICC DZ; red, ICC MZ > ICC DZ and p < 0.1 between ACE and CE models. EC1.4.1.9, EC2.6.1.42, EC4.2.1.10, EC1.1.1.282, EC1.1.1.25, EC5.4.99.5, EC4.2.1.51, EC2.6.1.1, EC2.6.1.57, and EC4.1.1.48 are bidirectional enzymes. EC4.2.1.9, EC4.3.1.19, EC4.2.1.19, EC3.1.3.15, and EC3.5.4.25 are mapped by two or more KOs for multiple metabolic reactions. The primary pathways for L-proline (EC2.7.2.11, EC1.2.1.41, and EC1.5.1.2) and L-arginine (EC6.3.4.5 and EC3.5.3.6) biosynthesis are significantly heritable (B). More detailed data are available in Table S4A for individual KOs.

Figure 5. Heritability of T2D MLGs and Butyrate Biosynthesis Pathways

(A) MLGs (>100 genes) from Qin et al. (2012) were profiled in the TwinsUK cohort. Blue nodes, present in less than 50% of the samples; green nodes, ICC MZ < DZ; pink nodes, ICC MZ > DZ; red nodes, ICC MZ > DZ and p < 0.1 between ACE and CE models (Table S5A). Light to dark-blue edges, Spearman’s correlation >0.4; light to dark-red edges, Spearman’s correlation < −0.4. (B) Pathways for butyrate biosynthesis were drawn according to Vital et al. (2014). Green arrows, ICC MZ < DZ; pink arrows, ICC MZ > DZ; red arrows, ICC MZ > DZ and p < 0.1 between ACE and CE models (Table S5B). The black line from crotonoyl-CoA to butyryl-CoA was not analyzed, because it was shared by all four pathways leading to Butyryl-CoA.

Figure 6. Gut Microbiome SNPs Detected in the TwinsUK Cohort

(A) SNP density in the 152 reference bacterial genomes with a cumulative coverage of at least 10x in the 250 samples. The bacterial genomes were ordered according to the cumulative coverage (green circles, Table S6). The coverage in each sample (black circles) was also plotted, with the maximum coverage among samples highlighted in beige. (B, E, and H) SNP similarity score within lean, overweight, and obese groups, calculated from all 152 reference genomes (B), for A. muciniphila (E) or D. longicatena only (H). p values from Wilcoxon rank-sum tests. 130 subjects in the lean group have a normal BMI (18.50–24.99) except for four underweight subjects (<18.5). (C, F, and I) SNP similarity score between the different BMI groups, calculated from all 152 reference genomes (C), for A. muciniphila (F) or D. longicatena only (I). p values from Wilcoxon rank-sum tests. (D, G, and J) SNP similarity score within and between geographic regions (Figure 2) calculated from all 152 reference genomes (D), for A. muciniphila (G) or D. longicatena only (J). p values from Wilcoxon rank-sum tests.

Figure 7. Greater Sharing of Gut Microbiome SNPs between Twins

(A) SNP similarity score between twins compared to unrelated samples, calculated from all 152 reference genomes. p = 0.0142 between MZ and unpaired samples, p = 0.1805 between MZ and DZ, p = 0.0328 between DZ and unpaired, one-tailed Wilcoxon rank-sum test. (B and C) SNP similarity score between MZ twins compared to DZ twins for A. shahii (B) and A. putredinis (C). p = 0.0408 and 0.01144, respectively, one-tailed Wilcoxon rank-sum test.