Large-scale comparative metagenomics of Blastocystis, a common member of the human gut microbiome

Abstract

The influence of unicellular eukaryotic microorganisms on human gut health and disease is still largely unexplored. Blastocystis spp. commonly colonize the gut, but its clinical significance and ecological role are currently unsettled. We have developed a high-sensitivity bioinformatic pipeline to detect Blastocystis subtypes (STs) from shotgun metagenomics, and applied it to 12 large data sets, comprising 1689 subjects of different geographic origin, disease status and lifestyle. We confirmed and extended previous observations on the high prevalence the microrganism in the population (14.9%), its non-random and ST-specific distribution, and its ability to cause persistent (asymptomatic) colonization. These findings, along with the higher prevalence observed in non-westernized individuals, the lack of positive association with any of the disease considered, and decreased presence in individuals with dysbiosis associated with colorectal cancer and Crohn's disease, strongly suggest that Blastocystis is a component of the healthy gut microbiome. Further, we found an inverse association between body mass index and Blastocystis, and strong co-occurrence with archaeal organisms (Methanobrevibacter smithii) and several bacterial species. The association of specific microbial community structures with Blastocystis was confirmed by the high predictability (up to 0.91 area under the curve) of the microorganism colonization based on the species-level composition of the microbiome. Finally, we reconstructed and functionally profiled 43 new draft Blastocystis genomes and discovered a higher intra subtype variability of ST1 and ST2 compared with ST3 and ST4. Altogether, we provide an in-depth epidemiologic, ecological, and genomic analysis of Blastocystis, and show how metagenomics can be crucial to advance population genomics of human parasites.

Figure 1

Prevalence of Blastocystis and Blastocystis subtypes in the different data sets, different continents, different European states, and between westernized and non-westernized subjects (see Supplementary Table 2 for more details). Stacked barplots show the prevalence in each category; numbers below the bars refer to the number of samples in the corresponding category, where duplicates from the same subject are eventually discarded. Statistical significance was assessed by Fisher's exact test.

Figure 2

Blastocystis prevalence in BMI classes (a) and different health conditions (b) for the considered data sets. Barplots show the prevalence of Blastocystis in different health conditions reported in the analyzed data sets. BMI classes considered were underweight (Un), normal (No), overweight (Ov) and obese (Ob). The total number of samples in each class and data set is reported below the bars. Bars associated with a total number of samples less than four are not shown. Note that scales in panels A and B are different. Abbreviations: CD, Crohn’s disease; IGT, impaired glucose tolerance; T2D, type 2 diabetes; UC, ulcerative colitis. Fisher's exact test was used as statistical significance test.

Figure 3

Breadth of coverage (a) and relative abundance (b) of Blastocystis in subjects colonized over two timepoints (see Supplementary Table 3 for more details). In the breadth of coverage plots, samples below the threshold of detection are also indicated. The breadth of coverage represents the fraction of the reference genome covered by at least one metagenomic read. The relative abundance is estimated by dividing the number of reads mapped to the Blastocystis reference genome with the total number of reads in the sample.

Figure 4

Phylogenetic relation between the 9 available Blastocystis reference genomes and 43 newly reconstructed genome assemblies from metagenomes. From the overall phylogenetic tree (a) we also report the subtrees of the four subtypes with more than 3 genomes (b–e) and compare the sequence diversity they span (f). Maximum likelihood phylogenetic trees were inferred using concatenated aligned shared genomic regions identified in reference genomes and assemblies (see Materials and methods). The asterisk highlights samples acquired from the same subject at two different timepoints. Black filled circles denote bootstrap support greater than 80%. The scale bar represents the average SNV rate calculated on the pairwise alignment.

Figure 5

Functional annotation analysis of the 19 reconstructed genomes spanning four Blastocystis STs. Less than half of the genes predicted by MAKER (Cantarel et al., 2008) were assigned to known COG functional categories using the eggNOG (Huerta-Cepas et al., 2016) database (a, see Materials and methods). These annotated genes can be grouped into 23 broad COG categories (b) that are a variable fraction of the total annotated genes. The asterisks denote categories for which one-way ANOVA statistical test gave P<1E−04. Hierarchical clustering performed on the more specific KOG functions show that samples associated with the same ST cluster together (c).

Figure 6

The presence (or absence) of Blastocystis is associated with major differences in the intestinal microbiome. Some species are strongly associated with the presence (a–c) or absence (d) of Blastocystis (plots for additional microbes are reported in Supplementary Figure 7). Boxplots report the distribution of abundances in samples with and without Blastocystis. Blue asterisks denote data sets where significant differences exist between the absence and presence of Blastocystis. LEfSe analysis (e) showed several other microorganisms statistically associated (α=0.05) with Blastocystis presence at high effect size (threshold at 3.3). Machine learning-based approach reveals that the microbiome signature is predictive for the presence of Blastocystis (f). This is valid not only when considering specific data sets with an unbiased cross-validation procedure, but also when predicting the presence of the parasite in a given data set considering only the samples from other independent studies (leave-one-data set-out approach).