Blastocystis Colonization Is Associated with Increased Diversity and Altered Gut Bacterial Communities in Healthy Malian Children

Blastocystis is the most common protozoan colonizing the gut of vertebrates. It modulates the human digestive microbiota in the absence of inflammation and gastrointestinal disease. Although it has been associated with human diseases, including inflammatory bowel disease, its pathogenicity remains controversial. This study aimed to assess the influence of Blastocystis on the gut bacterial communities in healthy children. We conducted a cross-sectional study on 147 Blastocystis-colonized and 149 Blastocystis-noncolonized Malian children, with Blastocystis colonization assessed by real-time PCR and gut microbial communities characterized via 16S rRNA gene (Illumina MiSeq) sequencing and bioinformatics analysis. The gut microbiota diversity was higher in Blastocystis-colonized compared to Blastocystis-noncolonized children. The phyla Firmicutes, Elusimicrobia, Lentisphaerae, and Euryarchaeota were higher in Blastocystis-colonized children, whereas Actinobacteria, Proteobacteria, unassigned bacteria, and Deinococcus–Thermus were higher in Blastocystis-noncolonized children. Moreover, Faecalibacterium prausnitzii (family Ruminococcaceae) and Roseburia sp. (family Lachnospiraceae) abundance was higher in Blastocystis-colonized children. We conclude that Blastocystis colonization is significantly associated with a higher diversity of the gut bacterial communities in healthy children, while it is not associated with the presence of potentially pathogenic bacteria in the human gut.


Introduction
Blastocystis is a genus of unicellular protozoan of the stramenopile group in the eukaryotic domain, described in 1911 by Alexeieff [1]. Many types of Blastocystis live in anaerobic conditions within the 24 h, solid stools were diluted v/v with 10X phosphate-buffered saline (PBS pH 7.4, RNase-free); all stools were aliquoted into 1 mL tubes and stored at −20 • C before being shipped to Marseille.

DNA Extraction
The total stool DNA was extracted with the EZ1 DNA tissue kit (Qiagen GmbH, Hilden, Germany) according to the procedure described herein [28]. Briefly, 200 mg (200 µL if liquid stools) and 350 µL of G2 lysis buffer were added in a 1.5 mL tube containing 200 mg of 2 mm glass powder. The sample was grinded with the FastPrep-24 TM 5G V. 6005.1 (M.P Biomedicals, LLC Santa Ana CA, USA) at 6 m/s for 40 s and then incubated at 100 • C for 10 min. After centrifugation at 10,000 × g for 10 min, 200 µL of the supernatant was recovered, and 10 µL of proteinase K was added followed by overnight incubation at 55 • C. We added 10 µL of a synthetic sequence of 142 bp (5 GCTACTGAGTCGTACCTAATGCATGACCTAGAGCACTCGACTGTTTATCAGTGTCG AGACTCGACGCATGCACGTACGAACCTAGCTGTCAGCAATCGCGAATGCCTACTAAGTAGCG AACTTTAGCGAATCGCGATACGAC-3 ) as a first extraction control at the concentration of 200 nmol diluted at 10 −10 in supernatants after the proteinase K digestion. The automated EZ1 Advanced XL (QIAGEN Instruments, Hombrechtikon, Switzerland) with the DNA card bacteria V 1.066069118 QIAGEN and the EZ1 DNA tissue kit was used to obtain 200 µL total DNA according to the procedure of extraction. Real-time PCR was performed on each DNA extracted for amplification of a synthetic sequence using the primers TissF_5 -CTGAGTCGTACCTAATGCATGACC-3 and TissR_5 -GTATCGCGATTCGCTAAAGTTC-3 and the probe TissP_6FAM-5 -TCGAGACTCGACGCATGCACG-Tamra-3 . The second extraction control of bacterial DNA was also carried out on the total DNA as described in [29]. The lysis buffer used for the extraction was used for the PCR negative controls. The extracted DNA was kept at 4 • C and immediately used as template for the PCR detection of gastrointestinal eukaryotic pathogens.

Real-Time PCR Assays
The primers and probes detailed in the paper of Sow et al. [30] were used to assay the twenty gastrointestinal eukaryotic pathogens on the CFX96TM and CFX384TM real-time PCR detection systems (BIO-RAD, Life Science, Marnes-la-coquette, France). The amplification reaction included 10 µL of master mix (Roche diagnostics GmbH, Mannheim, Germany), 0.5 µL of each primer, 0.5 µL of probe, 3 µL of distilled water, 0.5 µL of uracil-DNA glycosylase (UDG), and 5 µL of DNA in a total volume of 20 µL. The amplification program was as follows: 2 min at 50 • C, 5 min at 95 • C, followed by 40 cycles of 5 s at 95 • C and 30 s at 60 • C. A parasite-specific plasmid and a mixture of the DNA-free amplification reaction were used as positive and negative controls, respectively. Analyzed samples with no detectable amplification and Ct values greater than 39 were considered negative.

DNA Extraction and 16S Metabarcoding
DNA was extracted from stool samples after a first mechanical lysis step performed with acid-washed powder (≤106 µm) glass beads (G4649-500G Sigma-Aldrich, St. Quentin Fallavier, France) and 0.5 mm glass bead cell disruption media (Scientific Industries Inc., Bohemia, NY, USA) using a FastPrep BIO 101 instrument (Qbiogene, Strasbourg, France) at maximum speed (6.5 m/s) for 90s. Then, the stools were further lysed by two methods: (1) the classical lysis and protease step followed by purification on the NucleoSpin tissue kit (Macherey Nagel, Hoerdt, France) (protocol 1) and (2) deglycosylation and purification on the EZ1 Advanced XL device (Qiagen, Courtaboeuf, France) (protocol 5) [31]. Samples were first amplified on each of these two extraction products, then pooled and barcoded. The 16S rRNA sequencing was performed on the MiSeq system (Illumina, Inc, San Diego, CA, USA) with paired-end strategy, constructed according to the 16S metagenomic sequencing library preparation (Illumina). For each metagenomic DNA extraction protocol, the 16S "V3-V4" regions was amplified by PCR for 45 cycles using the Kapa HiFi Hotstart ReadyMix 2x (Kapa Biosystems Inc, Wilmington, MA, USA) and the surrounding conserved region V3_V4 primers with overhang adapters (FwOvAd_341F TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGCCTACGGGNGGCWGCAG; RevOvAd_785R GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGGACTACHVGGGTATCTAATCC). After amplicon purification on AMPure beads (Beckman Coulter Inc., Fullerton, CA, USA), DNA concentration was measured using high-sensitivity Qubit technology (Beckman Coulter Inc, Fullerton, CA, USA) and then diluted to 3.5 ng/ µL. At this step, the library of protocol 1 was pooled volume to volume with the library of protocol 5 so that 15 ng was involved in a subsequent limited cycle PCR, where Illumina sequencing adapters and dual-index barcodes were added to the amplicon. After purification on AMPure beads (Beckman Coulter Inc, Fullerton, CA, USA), this library was pooled with 95 other multiplexed samples. The global concentration was quantified by a Qubit assay with the high-sensitivity kit (Life technologies, Carlsbad, CA, USA). Before loading for sequencing on a MiSeq system (Illumina Inc., San Diego, CA, USA), the pool was diluted at 8 pM. Automated cluster generation and paired-end sequencing with dual index reads were performed in a single 39 h run in 2 × 250 bp. The paired reads were filtered according to the read qualities. The raw data were configured in FASTAQ files for R1 and R2 reads.

Bioinformatics Analysis
The 16S metabarcode bioinformatic analysis was performed with the MetaGX tool according to the bioinformatic company XEGEN's protocol [32]. Each sample read was analyzed by the VSEARCH tool to check raw data quality. The data cleaning was mainly done with the MiSeq sequencer software and then with the Cutadapt tool to eliminate sequencing primers and reads that were poor in quality and/or too short. The paired reads were kept and merged by the PANDAseq tool. The QIMME tool was used for merging and labeling the samples and then for clustering and operational taxonomic unit (OTU) creation using 97% of the clustering threshold. The QIMME tool allowed the OTUs to be filtered and a representative sequence for each OTU to be chosen.
The taxonomic assignment for each OTU was performed via BLAST querying the SILVA and the IHU MI culturomics in-house databases. Blast hits were sorted on the basis of coverage and percentage of identity depending on the (1) presence of one or more blast hits associated with a reference sequence (100% coverage; identity >97% corresponds to OTU's assignment to the species associated with the best blast hit); (2) presence of less relevant blast hits (identity between 95% and 97%: assignment to genus level; between 90% and 95%: assignment to the family; below 90%: assignment to the kingdom) with creation of a putative species in each case; (3) no blast hits (creation of putative new bacterial species).

Statistical Analysis
Multivariate linear regression analysis was used to test the association between Blastocystis and sociodemographic and clinical covariates. The number of observed OTUs; the Shannon, Simpson, and Chao-1 diversity indices; the abundance of OTUs; and the Bray-Curtis dissimilarity index were computed with the PAST3 software package (PAleontological Statistics Sofwware Version 3.20 [https://palaeo-electronica.org/2001_1/past/issue1_01.htm]). The median, interquartile range, mean, and standard deviation of these covariates were tabulated. Alpha diversity in Blastocystis-colonized and Blastocystis-noncolonized by were compared via the Mann-Whitney-Wilcoxon test with the GraphPad Prism software version 7.00 for Windows. Principal coordinate analysis (PCoA) in PAST3 was used to display the association of bacterial communities with Blastocystis-colonized or Blastocystis-noncolonized children. The bacterial community structure of the two groups was compared using the nonparametric statistical analysis of significant difference between two (permutational multivariate analysis of variance (PERMANOVA)) tests. A linear discriminant effect size (LEfSe) analysis was used to estimate taxon effect values between the two groups (http://huttenhower.sph.harvard.edu/lefse/). All statistical tests were two-sided, and p < 0.05 was considered statistically significant.

Ethical Considerations
The study was conducted in accordance with the Declaration of Helsinki, and the protocol was approved by the Ethics Committee of the Faculty of Medicine of Mali (No 2017/133/CE/FMPOS).
Written informed consent for participation in the study was obtained from each child and at least one of his/her parents or responsible person, and written assent was obtained from older children.

Characteristics of Study Subjects
In this study, we included 300 healthy participants from a prospective malaria cohort study with a mean age of 8 years, consisting of 154 (51.3%) females and 146 (48.7%) males (p = 0.97). We collected one stool sample from each child. Four samples were excluded, giving a total of 296 (98.7%) samples. Real-time PCR detected 147/296 (49.7%) Blastocystis-positive stool samples. The mean age of Blastocystis-positive and Blastocystis-negative subjects was 7.5 years and 8.3 years, respectively (p = 0.32). The male/female sex ratio was 0.97 and 0.89 in Blastocystis-positive and Blastocystis-negative subjects, respectively (p = 0.93).
Multivariate analysis highlighted a statistically significant association of Blastocystis colonization with age (p = 0.001) and weight (p = 0.03). No other sociodemographic, clinical, or biological characteristic was statistically significantly associated with Blastocystis colonization (Table S1).
There was a statistically significantly higher OTU richness and bacterial diversity in Blastocystis-colonized compared to Blastocystis-noncolonized children. Indeed, both Shannon and Simpson diversity indices showed that the diversity of the gut bacterial communities in Blastocystisnoncolonized children was lower than in Blastocystis-colonized subjects (p < 0.01) (Figure 1; Table 1). In contrast, both the Chao-1 and OTU richness indices showed that the richness of the gut bacterial communities were higher in Blastocystis-colonized children compared to Blastocystis-noncolonized children (p < 0.01) (Figure 1; Table 1). The Bray-Curtis dissimilarity index, which assessed the differences in bacterial community structure between Blastocystis-colonized and Blastocystis-noncolonized groups of children, was used in PCoA and PERMANOVA. The PCoA showed a clustering of the samples depending on the Blastocystis colonization status of the children; coordinate 1 (42.5%) and coordinate 2 (8.7%) scores explained 50% of the variance of the data (Figure 2). Furthermore, PERMANOVA showed a statistically significant difference in the bacterial community structure between the two groups (p = 0.0001).

Impact of Blastocystis on Gut Bacterial Communities
The heterogeneity of bacterial communities between Blastocystis-colonized and Blastocystisnoncolonized children was evaluated via LEfSe [33]. The results showed that the phyla Firmicutes, Elusimicrobia, Lentisphaerae, Euryarchaeota, and IHU_PP_Bacteria were significantly more abundant in the Blastocystis-colonized group than Actinobacteria, Proteobacteria, unassigned bacteria, and Deinococcus-Thermus, which were overrepresented in the Blastocystis-noncolonized children (Supplementary file 2: Figure S5A).

Impact of Blastocystis on Gut Bacterial Communities
The heterogeneity of bacterial communities between Blastocystis-colonized and Blastocystisnoncolonized children was evaluated via LEfSe [33]. The results showed that the phyla Firmicutes, Elusimicrobia, Lentisphaerae, Euryarchaeota, and IHU_PP_Bacteria were significantly more abundant in the Blastocystis-colonized group than Actinobacteria, Proteobacteria, unassigned

Impact of Blastocystis on Gut Bacterial Communities
The heterogeneity of bacterial communities between Blastocystis-colonized and Blastocystisnoncolonized children was evaluated via LEfSe [33]. The results showed that the phyla Firmicutes, Elusimicrobia, Lentisphaerae, Euryarchaeota, and IHU_PP_Bacteria were significantly more abundant in the Blastocystis-colonized group than Actinobacteria, Proteobacteria, unassigned bacteria, and Deinococcus-Thermus, which were overrepresented in the Blastocystis-noncolonized children (Supplementary file 2: Figure S5A).
At the genus level, Ruminococcus and Clostridium were among the most abundant genera in the Blastocystis-colonized children, whereas Streptococcus, Bifidobacterium, and Shigella were more abundant in Blastocystis-noncolonized children (Supplementary file 2: Figure S5E).
At the species level, Clostridium saudii, Methanobrevibacter smithii, and a few other species were the most abundant in the Blastocystis-colonized children, whereas more Streptococcus sp., Bifidobacterium sp., Shigella sp., and a few other species were the most abundant in Blastocystis-noncolonized children (p < 0.05). Some species were notably abundant in the group of Blastocystis-noncolonized children compared to Blastocystis-colonized children (p < 0.05). (Supplementary file2: Figure S5F).

Discussion
This study's main finding was that both diversity and richness of gut bacterial communities was higher in healthy Blastocystis-colonized children compared to Blastocystis-noncolonized ones. A limitation of our study was that the influence of the children's lifestyle and diet on the bacterial communities was not assessed. The major strength of our study, compared with previous ones [34], was the robustness of our findings, which was supported by the relatively large sample size of 296 healthy Malian children.
It has been reported that the diversity of the eukaryotic gut communities in healthy humans is of relatively low abundance, stable over time, and dominated by Blastocystis subtypes [17]. The influence of Blastocystis on the gut bacterial communities has been assessed in several studies [34][35][36][37][38][39][40]. In line with our findings, many studies [34,39,41] have found that Blastocystis colonization is associated with a relatively increased diversity of gut bacterial communities. In particular, the diversity was relatively higher among Swedish travelers colonized with Blastocystis, which was also associated with both a "healthy" gut microbiota profile and a vegetable-rich diet [37]. Furthermore, a high bacterial community diversity has been observed in Blastocystis-colonized patients [34]. A study evaluating the influence of Giardia duodenalis, Entamoeba spp., and Blastocystis hominis infections on the structure of bacterial communities in symptomatic and asymptomatic subjects (adults/infants) in Côte d'Ivoire found that the Faecalibacterium prausnitzii/Escherichia coli ratio increased in the subjects carrying Entamoeba spp. and Blastocystis hominis compared to those with no enteric protozoan [39]. In contrast, the diversity of bacterial community decreased in Blastocystis-colonized patients with hepatic encephalopathy [35]. However, a decreased bacterial diversity has been observed in metabolic or infectious diseases, such as IBD or infection with enteric pathogens, which are associated with inflammation of the lower gastrointestinal tract [21,22,42].
Our study showed that the phyla Firmicutes and Bacteroidetes were most abundant followed by Actinobacteria and Proteobacteria in all children. This trend was confirmed in a study investigating the association of Blastocystis carried by patients with irritable bowel syndrome (IBS) with bacterial communities [38]. The differences in the bacterial community structure of Blastocystis-colonized and Blastocystis-noncolonized children explained 50% of the variance in a coordinated main component analysis. Furthermore, there was a statistically significant difference (p = 0.0001) using a nonparametric similarity test (PERMANOVA) with the Bray-Curtis similarity measure, supporting the disparity of the bacterial community structures in Blastocystis-colonized and Blastocystis-noncolonized children. Similarly, Audebert et al. found that Blastocystis colonization explained 30% of the variance in bacterial community structure, which was statistically significant (p = 0.0001) via PERMANOVA [34].
High-throughput Illumina MiSeq sequencing technology was used to characterize the gut bacterial community structure of the Blastocystis-colonized and Blastocystis-noncolonized children. The LEfSe, designed to discriminate the taxa associated with each group, highlighted a significant increase in the phyla Firmicutes, Elusimicrobia, Lentisphaerae, and Euryarchaeota in Blastocystis-colonized children, while Actinobacteria, Proteobacteria, unassigned bacteria, and Deinococcus-Thermus were higher in Blastocystis-noncolonized children. Therefore, Blastocystis colonization was associated with eubiosis, a condition characterized by a higher proportion of "beneficial bacteria" (Firmicutes and Bacteroides) than "probable pathogenic bacteria" (Proteobacteria). The authors hypothesized that Blastocystis colonization could be a surrogate marker of a "healthy gut microbiota". In line with this proposition, it has been reported that Blastocystis colonization is more frequent in healthy subjects compared to patients with IBD [18,19]. Furthermore, we found higher bacterial species diversity in Blastocystis-colonized compared to Blastocystis-noncolonized children. LEfSe highlighted that Clostridiaceae, Ruminococcaceae, and Lachnospiraceae were more abundant in Blastocystis-colonized children, whereas Streptococcaceae, Bifidobacteriaceae, Enterobacteriaceae, and Leuconostocaceae were more abundant in Blastocystis-noncolonized children. More precisely, Faecalibacterium prausnitzii (family Ruminococcaceae) and Roseburia sp. (family Lachnospiraceae) were relatively more abundant in children colonized by Blastocystis. This is of particular interest because other authors have reported that Faecalibacterium prausnitzii and Roseburia sp. abundance is decreased in the gut bacterial community of Crohn's disease patients [43]. Both Faecalibacterium prausnitzii and Roseburia sp. play an important protective role in gut physiology. They digest dietary fiber into short-chain fatty acids, especially butyrate, which provides an energy source for intestinal cells [44]. They also possess anti-inflammatory properties [45][46][47]. Similar to our findings, another study [34] showed that Faecalibacterium prausnitzii and Roseburia sp. were relatively more abundant in patients colonized by Blastocystis; this suggests that Blastocystis colonization could be used as a surrogate marker of a healthy gut microbiota.

Conclusions
Blastocystis colonization was significantly associated with a higher diversity and richness of the gut bacterial communities in healthy children. Also, Blastocystis colonization was associated with a higher proportion of "beneficial bacteria" (Firmicutes and Bacteroides) than "probable pathogenic bacteria" (Proteobacteria) in the human gut. Indeed, because of its hypothetical capacity to promote, or at least be associated with, an increased diversity of gut bacterial communities, documenting the presence of Blastocystis within a healthy microbiota donor would be critical in the management of patients whose diseases require fecal transplantation. Further studies aimed at elucidating the mechanisms by which Blastocystis influences the gut bacterial communities are warranted.
Author Contributions: S.R. and A.K. designed and conceived the study. D.C., A.K.K., S.K., S.D., and M.A.T. included the study participants and performed clinical and biological evaluations, sample collection, and data management. A.K., A.G., S.R., and F.B. did the statistical analysis. The manuscript was drafted by A.K. and edited by F.G., D.R., and S.R. All authors read and approved the final version of the manuscript.
Funding: This study was supported by funding from the IHU-Mediterranean Infection Foundation and the MARCAD DELTAS Africa Initiative grant DEL-15-10. The DELTAS Africa Initiative is an independent funding scheme of the African Academy of Sciences' (AAS) Alliance for Accelerating Excellence in Science in Africa (AESA) and supported by the New Partnership for Africa's Development Planning and Coordinating Agency (NEPAD Agency) with funding from the Wellcome Trust grant 107741/A/15/Z and the UK government. The views expressed in this publication are those of the author(s) and not necessarily those of AAS, NEPAD Agency, Wellcome Trust, or the UK government.