Design and application of a novel two-amplicon approach for defining eukaryotic microbiota

Abstract

Background: Due to a lack of systematic diagnostics, our understanding of the diversity and role of eukaryotic microbiota in human health is limited. While studies have shown fungal communities to be significant modulators of human health, information on the prevalence of taxa such as protozoa and helminths has been limited to a small number of species for which targeted molecular diagnostics are available. To probe the diversity of eukaryotic microbes and their relationships with other members of the microbiota, we applied in silico and experimental approaches to design a novel two-amplicon surveillance tool, based on sequencing regions of ribosomal RNA genes and their internal transcribed spacers. We subsequently demonstrated the utility of our approach by characterizing the eukaryotic microbiota of 46 hospitalized Malawian children suffering from Severe Acute Malnutrition (SAM).

Results: Through in silico analysis and validation on a diverse panel of eukaryotes, we identified 18S rRNA variable genetic regions 4 and 5 (18S V4 V5), together with a region encoding 28S rRNA variable genetic region 2 and the internal transcribed spacers (transITS), as optimal for the systematic classification of eukaryotes. Sequencing of these regions revealed protozoa in all stool samples from children with SAM and helminths in most, including several eukaryotes previously implicated in malnutrition and diarrheal disease. Clinical comparisons revealed no association between protozoan parasites and diarrhea or HIV reactivity. However, the presence of Blastocystis correlated with bacterial alpha diversity and increased abundance of specific taxa, including Sporobacter, Cellulosibacter, Oscillibacter, and Roseburia.

Conclusion: We suggest this novel two-amplicon based strategy will prove an effective tool to deliver new insights into the role of eukaryotic microbiota in health and disease.

Keywords: Amplicon sequencing; Eukaryotes; Malnutrition; Microbiome; Next-generation sequencing; Parasites.

Fig. 1

Identification and evaluation of biomarker regions in 18S and 28S rRNA genes. a Sequence variability of protozoan 18S and 28S rRNA genes. Variability was calculated as the Shannon entropy at each position along aligned protozoan sequences. Identified variable regions are numbered and those selected for further analysis are highlighted in purple. The graphic at the top depicts the canonical arrangement of eukaryotic rRNA genes. b Accuracy of taxonomic classification of variable regions to the phylum, genus, and species levels. Taxonomies of 100 copies of variable regions in each gene, randomly mutated with 1% frequency, were compared to the best match identified by BWA-MEM in the SILVA v128 database. Results are shown for protozoan (19,427 18S; 1061 28S), helminth (3882 18S; 1237 28S), and fungal (14,657 18S; 3625 28S) rRNA genes

Fig. 2

In silico and experimental evaluation of amplicon primers. Predicted DNA primer binding a to eukaryote, archaeal, and bacterial genes and b to protozoan and helminth kingdoms. Primers used in previous published 18S rRNA gene sequencing studies and the highest scoring primers designed for this study are shown for three potential amplicon regions. Tests were performed using SILVA TestProbe3.0 against the v128 Ref_Nr99 (18S) and Ref (28S) databases. c PCR amplification of the V4 V5 amplicon from the 18S rRNA gene using V4-1_F and V4-4_R, or 515f and 1119r primers (top), and the transITS amplicon using V9-9_F and V2-rev6_R primers (bottom). Reactions include a no template control, and template genomes from (1–15) T. gondii RH, N. caninum [NcBahia], B. hominis, T. vaginalis SD7, T. muris G1, T. cruzi [DA], L. panamensis [WR470], G. intestinalis [G3M], E. dispar, S. cerevisiae, C. albicans, A. fumigatus, C. elegans, S. mediterranea, and a human foreskin fibroblast cell line (ATCC SCRC-1041)

Fig. 3

Eukaryotic and bacterial microbiota of Malawian cohort suffering from SAM. a Relative sequence abundance of microbial eukaryotic phyla, identified by the 18S rRNA gene V4 V5 amplicon (n = 44). Columns represent individual patients. b Numbers of protozoan genera identified in patients from 18S rRNA gene V4 V5 data. c Bayesian phylogenetic trees showing relationships between sequence representatives from clinical samples and reference species. Posterior probability values indicate branch support. d Relative sequence abundance of microbial eukaryotic phyla, identified by transITS (n = 46) amplicons. e Similarity of taxonomic profiles generated by the two amplicon regions, quantified by co-occurrence of eukaryotic genera within patients compared to between patients. Shown are distributions of within-patient and between-patient normalized Hamming distances, averaged over a range of presence/absence thresholds (2–100 reads). f Bray-Curtis compositional dissimilarity of eukaryotic genera identified by the two amplicon regions within and between patients. Significance was tested using the Wilcoxon rank sum test. g Venn diagram showing the overlap of eukaryotic genera identified using the two amplicon regions. A subset of organisms of interest is indicated. h Relative sequence abundance of bacterial phyla identified in 16S rRNA gene amplicon data

Fig. 4

Classification profiles of the transITS amplicon. a Graphs show numbers of sequences (y-axis) classified to indicated eukaryotic taxa, at various minimum sequence coverage thresholds (x-axis). Colors indicate the maximum percent mismatch between amplicon and reference sequences. b Read coverage depicted across positions of the amplicon. Shown are total numbers of reads (violet) and reads with ≥ 97% sequence identity across ≥ 70% of the read length (blue) classified to a genus. Graphics below plots indicate available (red) or missing (gray) reference sequences for the genus

Fig. 5

Associations between eukaryotic and bacterial microbiota. a Co-occurrence analysis between eukaryotic and bacterial microbiota. Shown is a heatmap of Hamming distances between genera represented by ≥ 5 reads in 20–80% of samples (n = 44). Red indicates a high co-occurrence score. The detail highlights bacterial genera clustering with Blastocystis, and the yellow line denotes a cluster of Enterobacteriaceae. b Association of bacterial composition with the presence of Blastocystis. Plots were generated with ALDEx2 and show OTU fold change versus median abundance (left) and OTU fold change between versus within conditions (right). Red points indicate OTUs with significant compositional changes (Wilcoxon rank test). c sPLS-DA analysis showing separation of patient samples based on the abundance of bacterial genera in Blastocystis positive and negative patients. d Mean relative sequence abundances of a subset of bacterial genera, in Blastocystis positive and negative patients. Error bars represent standard error. e Comparison of alpha diversity of samples positive and negative for Blastocystis (15; 29), Cryptosporidium (12; 32), and Giardia (14; 26), measured using the Shannon index. The presence of Blastocystis and Cryptosporidium were counted with a minimum of five reads. *p < 0.05, **p < 0.005