Culturomics and Amplicon-based Metagenomic Approaches for the Study of Fungal Population in Human Gut Microbiota.

Herein, the mycobiota was characterized in fecal samples from sick patients and healthy subjects, collected from different geographical locations and using both culturomics and amplicon-based metagenomics approaches. Using the culturomics approach, a total of 17,800 fungal colonies were isolated from 14 fecal samples, and resulted in the isolation of 41 fungal species, of which 10 species had not been previously reported in the human gut. Deep sequencing of fungal-directed ITS1 and ITS2 amplicons led to the detection of a total of 142 OTUs and 173 OTUs from the ITS1 and ITS2 regions, respectively. Ascomycota composed the largest fraction of the total OTUs analyzed (78.9% and 68.2% of the OTUs from the ITS1 and ITS2 regions, respectively), followed by Basidiomycota (16.9% and 30.1% of the OTUs from the ITS1 and ITS2 regions, respectively). Interestingly, the results demonstrate that the ITS1/ITS2 amplicon sequencing provides different information about gut fungal communities compared to culturomics, though both approaches complete each other in assessing fungal diversity in fecal samples. We also report higher fungal diversity and abundance in patients compared to healthy subjects. In conclusion, combining both culturomic and amplicon-based metagenomic approaches may be a novel strategy towards analyzing fungal compositions in the human gut.

Introduction

The gastrointestinal tract is a complex human organ that is exposed to the external environment, and which acts as a semipermeable multilayer system for the absorption of the nutrients and macromolecules required for human metabolic processes[1]. The microbiota has been described in the literature as playing essential roles in metabolism, immune function and gene expression[2]. Moreover, the microbiota’s potential involvement in diseases such as allergies, inflammatory bowel disease, and metabolic and degenerative disorders has also been described[3]. Human gut microbiota is composed of bacteria, Archaea, and small eukaryotes[4,5]. Bacteria represent the majority of these microbial communities inhabiting the human gut, with four major phyla resident (namely, Firmicutes, Bacteroidetes, Actinobacteria and Proteobacteria)[6]. Unlike the human gut’s bacterial community, the composition and the diversity of eukaryotes remains largely unexplored due to their relative low abundance in the human gut or their neglect by both culture and molecular analyses[7]. Fungi represent an important component of the eukaryotic microorganisms in the human gastrointestinal tract[8,9]. Although the presence of fungal species in the human gut appears to be mainly commensal or mutualistic[10], some yeast and filamentous fungi are responsible for emerging pathologies, particularly in immunocompromised patients[11]. Fungal species detected in the human body belong mainly to three different phyla: Ascomycota, Basidomycota and Zygomycota[12]. Several techniques are used for fungal identification, including culture, microscopy, sequencing and, more recently, Matrix Assisted Laser Desorption/ionization Time-Of-Flight (MALDI-TOF) mass spectrometry[13] and metagenomics[9,12,14-20]. A total of 278 fungal species have been reported in the human gut[13]. These species have been identified either by culture-dependent or culture-independent methods[13]. Among these fungal species, 179/278 (64%) are classified in the Ascomycota phylum, 83/278 (30%) in the Basidiomycota phylum, with 7/278 (3%) belonging to Zygomycota and 5/278 (2%) species in the Microsporidia group. However, 4/278 (1%) of species are still unclassified. Recent molecular studies have demonstrated the efficiency of these techniques in identifying a high number of fungi[13]. Notable progress has been made in gut fungal diversity with the advent of next generation sequencing (NGS)[21], which has led to the expansion of the knowledge of digestive tract mycobiota and has enhanced the understanding of the roles of these fungi in healthy states and in diseases[15-20,22]. Culture-dependent methods allow the identification of only 75 (27%) fungal species in human gut microbiota, using mainly fecal samples[13]. Culturomics has emerged recently as a successful tool to isolate high number of bacteria and to identify new species[23-25]. It has shifted the view on human gut microbiota by supplying “new” bacterial diversity not previously captured by metagenomics. However, culturomic methodology has not yet been comprehensively applied to describe the fungal population in the human gut. Thus, the objective of this study is to shed some light on fungal diversity from stools sampled from different geographic locations and from sick patients and healthy subjects using both culturomics and ITS1/ITS2 amplicon sequencing approaches.

Results

Fungal identification and species diversity

Positive cultures were obtained from all 14 human fecal samples used in this study. In a total, 17,800 fungal colonies were isolated and analyzed using MALDI-TOF mass spectrometry identification and ITS sequencing. Forty-one fungal species were identified (Supplementary Table 1).

Among the 41 fungal species isolated by culturomics, 16 (39%) could be identified successfully by MALDI-TOF mass spectrometry with spectra score ≥1.9 (Supplementary Table 1). All of the 16 identified fungal species were yeast strains. Saccharomyces cerevisiae was identified with the highest score (2.4). Many fungi, including mainly filamentous fungi (No = 25; 61%), could not be identified by mass spectrometry, but only by direct ITS sequencing (Supplementary Table 1).

Moreover, all 41 (100%) fungal strains were sequenced by targeting the ITS region of the ribosomal RNA gene (Supplementary Table 1). The ITS results confirmed those found by mass spectrometry for the yeast strains.

The fungal species isolated in this study belonged to 3 different phyla, with the majority (34/41; 82.9%) belonging to the Ascomycota phylum, 5/41 (12.2%) belonging to the Basidiomycota phylum and the remaining 2/41 (4.9%) of fungal species belonging to Zygomycota.

Fungal isolation according to culture media

A total of five culture media (Supplementary Table 2) and 12 conditions were used in this study. All strains were isolated under aerobic conditions, at 28 °C and at room temperature (22 °C). The best medium was Dixon agar, allowing the isolation of 30/41 (73.2%) strains, followed by modified Schadler Agar allowing the isolation of 22/41 (53.7%) of strains. Sabouraud medium allowed only the isolation of eight species (19.5%). Finally, BAM and PDA allowed the isolation of 6/41 (14.6%) and 5/41 (12.2%) fungi respectively (Supplementary Figure 1).

A pre-incubation step in liquid medium (containing sterile rumen fluid and sheep’s blood) was also used in this study. Although this step did not allow the isolation of more strains than those found directly in solid axenic media without pre-incubation, it resulted in the reduction of time spent to isolate strains, especially for filamentous species. Further, it let to more colonies for each fungus cultured on solid agar media.

Fungal diversity detected by Culturomics

Among the 14 fecal samples used in this study, half were collected from sick individuals including obese, HIV and marasmus patients. The number of fungal species per sample ranged between 2 to 10 and 3 to 12 in the healthy and sick patient group respectively. The most common isolated fungi were Aspergillus niger and Candida glabrata (both were isolated in 7 samples) followed by Candida parapsilosis (6 samples), Clavispora lusitaniae and Debaryomyces hansenii (each were isolated in 5 samples) and Candida albicans (4 samples) (Supplementary Figure 2).

Composition of fungal diversity using ITS1 and ITS2 dataset

The fungal diversity of the14 stool samples collected from different geographic locations was estimated using high throughput sequencing technology. Over 1,326,225 paired-reads were collected from both ITS1 and ITS2 regions using MiSeq sequencing system (Illumina). Of these, 904,412 reads with high quality sequences were then selected for downstream analysis. Clustering of these reads using 97% of similarity resulted in a total of 142 OTUs and 173 OTUs for the ITS1 (Fig. 1 and Supplementary Figure 3) and ITS2 (Fig. 2 and Supplementary Figure 4) regions, respectively. Ascomycota composed the largest fraction of the total OTUs (78.9%) analyzed from the ITS1 region, followed by Basidiomycota (16.9%). Finally, 4.2% of the OTUs were ascribed to unclassified fungal phyla. Approximately 10.6%, 69%, 7%, and 2.8% of the OTUs were identified at the species, genus, family and order, respectively. The remaining OTUs (10.6%) were assigned to the level of class or higher (Fig. 1, Supplementary Figure 3 and Supplementary Table 3). Notably, two OTUs from genus Torulaspora (OTU30 and 32) were detected in all 14 fecal samples (Supplementary Table 3). Whereas C. albicans (OTU16) and OTUs belonging to genus Debaryomyces (OTU47 and OUT51) and genus Torulaspora (OTU31, OTU37 and OTU46) were identified in most (93%) of the fecal samples (Supplementary Table 3). The ascomycete yeast from genus Aureobasidium (OTU142) and the mold from genus Aspergillus (OTU59) and Penicillium (OTU95) were amplified in more than 70% of all fecal samples. Similar to the ITS1 dataset, members of Ascomycota constituted the majority (68.2%) of the total fungal OTUs detected in the ITS2 dataset followed by Basidiomycota (30%). Whereas Mucoromycota represented only a small fraction (0.6%) of the total OTUs analyzed from the ITS2 region (Fig. 2, Supplementary Figure 4 and Supplementary Table 4). Finally, 1.2% of the OTUs were recognized as unclassified fungal phyla. The results from analyzing the ITS2 dataset revealed that 11%, 55.5%, 9.2% and 12.2% of the total obtained OTUs were assigned at the level of species, genus, family and order, respectively. The remaining OTUs (12.1%) were identified in the ITS2 data set at the level of class or higher (Fig. 2, Supplementary Figure 4 and Supplementary Table 4). Species such as C. albicans (OTU18), Kluyveromyces marxianus (OTU35), S. cerevisiae (OTU36), Malassezia globosa (OTU85) and Wallemia muriae (OTU171) were detected in all fecal samples (Supplementary Table 4). OTUs belonging to genus Aspergillus (OTU39, OTU87), Aureobasidium (OTU173), Beauveria (OTU144), Cladosporium (OTU124), Debaryomyces (OTU21), Meyerozyma (OTU15, OTU127), Penicillium (OTU88), Malassezia (OTU128) and Trichosporon (OTU129) were also represented in all fecal samples. Further, each of Aspergillus (OTU40), Penicillium (OTU42), Trichosporon (OTU64) and Aureobasidium (OTU38) were identified in more than (80%) of all fecal samples (Supplementary Table 4).

Figure 1

Fungal OTUs (at > 0.1%) detected from the amplification of the ITS1 region in the fecal samples of healthy and unhealthy subjects.

Figure 2

Fungal OTUs (at > 0.1%) detected from the amplification of the ITS2 region in the fecal samples of healthy and unhealthy subjects.

Compared fungal richness and diversity among fecal samples

Rarefaction curves (distance levels of 0.03) for fungal community in all fecal samples initiate to reach asymptotic levels except for samples (HIV20 and Bedouin7) (Supplementary Figure 5). Thus, the sequencing capabilities were almost representing the number of different fungal taxa amplified in ITS1 dataset. As with rarefaction curves of the ITS1 dataset, the rarefaction curves for the ITS2 dataset also reached saturation (Supplementary Figure 6), indicating sufficient depth of sequencing of the different fungal communities amplified within the ITS2 dataset.

The indices of diversity and abundance of intestinal fungi in fecal samples using the ITS1 dataset are well illustrated in Supplementary Table 5. We used the Chao-1 estimator and Shannon diversity indices to estimate the diversity of fungal communities in the different fecal samples (Supplementary Table 5). The highest Chao-1 value for the fungal community was obtained in the Obese3 fecal sample, whereas the lowest value detected in the Bedouin7 fecal sample (Supplementary Table 5). Furthermore, Shannon diversity index values showed relatively high fungal diversity in the HIV5 fecal sample compared with the rest of the fecal samples (Supplementary Table 5). In the ITS2 dataset, the highest Chao-1 value for the intestinal fungal community was obtained in the HIV5 fecal sample and the lowest value detected in the Touareg4 fecal sample. Whereas Shannon diversity index values showed that the fungal diversity in the Obese3 fecal sample was relatively greater than the rest of the fecal samples (Supplementary Table 5). Overall, the alpha diversity, including Chao-1 richness estimation and Shannon’s diversity index, revealed significantly higher fungal diversity and abundance in sick patients compared to healthy subjects (Supplementary Figure 7).

Discussion

The objective behind this work was to evaluate the diversity of fungi in fecal samples from healthy subjects and sick patients, collected from different geographical locations and using both culturomic and ITS1/ITS2 amplicon-based metagenomic approaches[23]. Culturomic technique has never been applied extensively to fungi, although preliminary studies have been published[8,9,14].

Using this approach, a total 41 fungal species were isolated, 10 of which had never been found in the human gut. Further, 4 such fungal species (namely Sporormiella minima, Pichia mandshurica, Penicillium glandicola and Ascophaera apis) had not been previously isolated from the human body. Dixon agar (DIX), Sabouraud (SAB) and Potato Dextrose agar (PDA) are commonly used for fungal species isolation[14,26,27]. More than 73.2% (30 /41) of fungal species were isolated with DIX. In fact, many studies have demonstrated that DIX is a suitable culture media for fungal isolation and allows the cultivation of lipophilic fungal species such as species belonging to the genus Malassezia [8,28]. Few species have been isolated with SAB, Banana Agar Medium (BAM) and PDA culture media. Finally, 22/41 (53.7%) fungal species were isolated with the modified Schaedler Agar (MSA). As such, it may be concluded that MSA is a good culture media for fungal species as it allows the isolation of many fungal strains. Liquid enrichment through the addition of compounds such as rumen fluid and sheep’s blood[23] can also significantly reduce the time required to culture strains, as was demonstrated by the results of our study.

Even though the use of selective and/or enrichment culture conditions coupled with MALDI-TOF MS identification is a common technique for the identification of microorganisms in many laboratories[27,29-33], there are nevertheless major limitations in the identification of certain isolates by MALDI-TOF MS. This is believed to be due to either the insufficiency or absence of reference spectra library in such databases[30,32,34,35], as well as limitations associated with classic culture techniques for fungal isolation[36].

The application of NGS-based metagenomic sequencing to assess the fungal diversity in human gut is becoming more common[15-20,22]. The ITS region has been used as a default marker for most gut fungal metagenomic studies[17-20,22,37] for the last decade. Due to the short length reads generated by most high-throughput sequencing platforms, researchers have been forced to select either the ITS1 or the ITS2[17-20,22,37] sublocus in their fungal diversity analyses. However, studies have shown that targeting either the ITS1 or ITS2 region separately introduces fungal taxonomic biases[38,39]. Therefore, focusing on both the ITS1 and ITS2 subloci may provide a better taxonomic resolution on the assessment of fungal diversity in environmental samples[40-42]. As such, in the current study, the fungal gut components were studied in feces sampled from sick patients and healthy subjects inhabiting different geographical locations using deep sequencing of fungal-directed ITS1 and ITS2 amplicons. A total of 142 OTUs and 173 OTUs were detected in each of the ITS1 and ITS2 datasets respectively. Bringing together the results from analyzing both the ITS1 and ITS2 datasets, a total of 163 fungal taxa (each have Unique Species Hypotheses (SH) Accession number) were obtained, among which 52 were unique to the ITS1 dataset, 81 were unique to the ITS2 dataset, and only 30 fungal taxa were common to both (Supplementary Figure 8). It is clear from our study that the results from the ITS1 and ITS2 dataset complemented each other in assessing fungal diversity in fecal samples.

On the other hand, combining the results of both culturomic and amplicon metagenomic approaches, resulted in the identification of 181 fungal taxa, among which 18 taxa were specific to the culturomics dataset, 140 taxa were specific to the ITS1 and ITS2 dataset and only 23 were common to both (Fig. 3). The culturomic approach in this study provided different information about gut fungal diversity, compared to the amplicon metagenomic approach based on the ITS1/ITS2 region of fungal rDNA, though the two approaches were complementary to each other. The incapability of ITS1/ITS2 amplicon sequencing method to amplify all cultured fungi could be explained by the biases due to DNA extraction and PCR amplification in amplicon metagenomic method[43]. Moreover, the short Illumina reads (~200 bp) that were generated in both ITS1 and ITS2 datasets could not provide adequate information for fungal taxon assignment compared to complete ITS sequences which confer taxonomic assignment up to species level.

Figure 3

Total unique fungal OTUs obtained by compiling the OTUs from the amplicon metagenomic approach (ITS1 and ITS2 datasets) and from fungal species obtained by the culturomic approach. Green species indicate that fungi were identified only in healthy fecal samples in this study. Red species indicate that fungi were detected only in non-healthy fecal samples. Black species indicate that fungi were detected in both healthy and non-healthy fecal samples. The blue circle indicates that fungal taxa were recovered from ITS1, ITS2 and culturomics whereas the green circle indicates that fungal taxa were recovered from both ITS1 and ITS2. Finally, red and orange circles indicate that fungal taxa were recovered from ITS1 and ITS2 respectively.

At the level of gut fungal richness, this study shows the presence of significantly higher fungal richness in sick patients than healthy subjects using ITS1/ITS2 amplicon sequencing method (Supplementary Figure 7). These results are consistent with previously reported data on fungal diversity in patients with inflammatory bowel disease or chronic hepatitis B[44-46].

Many yeast, such as members of Candida, Debaryomyces, Galactomyces and Rhodotorula, and mold, such as Aspergillus and Penicillium, were isolated from various fecal samples in our study. Candida spp. are considered typical commensal fungal species of the gastrointestinal tract whereas Aspergillus spp. are ubiquitous in nature[47]. Cryptococcus spp. are endemic mycoses that also have worldwide distribution. Recent studies have shown that genus Candida and Saccharomyces are isolated frequently from clinical samples[37], and that Candida spp. are the predominant commensal fungal species, while Rhodotorula spp. and Cryptococcus spp. are considered transient species or exogenous pathogenic fungi[48].

In conclusion, we provide novel insights into the fungal composition of the human gut, using both culturomic and ITS1/ITS2 amplicon sequencing approaches. We have highlighted the complementarity of these methods in providing a clearer picture of the true fungal diversity in the human gastrointestinal tract. As such, we suggest that both approaches should be analyzed in parallel in any gut fungal diversity analysis.

Materials and Methods

The samples

The study sample consisted of 14 stool samples, collected from various geographic locations and stored directly in sterile tubes at –80 °C until further use. These samples include 2 fecal samples from healthy subjects living in Marseille, France (MS1 and SO1), 1 sample from each of the Amazonia (Amazonia1) and Tuareg populations (Touareg4), 3 stool samples from a village in the Southern Region of Saudi Arabia (Bedouin2, Bedouin4 and Bedouin7), 3 stool samples from each of obese (Obese1, Obese3 and Obese5) and from patients with HIV (HIV5, HIV16 and HIV20), and a single stool from patient with marasmus (Marasmus10) (Supplementary Table 6). No antibiotic, antifungal or anti-parasitic therapy had been administrated to any enrolled subject at the time of the sample collection or in the previous 2 months. In this study, the samples obtained from patients (i.e. marasmus, obese and HIV patients) were grouped together as non-healthy group in order to compare their fungal alpha diversity to the healthy group (the remaining samples) given the fact that patients with these three diseases are known to have reduced gut bacterial diversity[49-51]. However, due to the heterogeneity of the studied samples (samples obtained from different diseases and from different geographic locations), the fungal beta diversity was not analyzed. All experiments and methods were performed in accordance with relevant guidelines and regulations. Written informed consent was obtained from all participants. The study was approved by the local ethics committee of the Institute Hospitalo-Universitare Méditerranée Infection Marseille, France (agreement number 2016–011).

Culturomics

In this study, we used culturomic[23] procedures developed previously by Lagier et al. for bacteria. Briefly, one g. of each fecal sample was diluted with 900 µl of Dulbecco’s Phosphate-Buffered Saline (DPBS) for solid and liquid fungal culture. 1/10th and 1/100th dilutions were performed. A total volume of 50 µl of each dilution was used for culturing on 5 solid culture media supplemented with 3 antibiotics; namely, colistin (30 mg/l), vancomycin (30 mg/l) and imipenem (30 mg/l). Solid culture media included Sabouraud agar (SAB)(Sigma-Aldrich, Saint-Quentin Fallavier, France), Dixon agar (DIX)[8,14,26], Potatoes Dextrose agar (PDA) (Sigma-Aldrich), modified Schaedler agar (MSA)(Sigma-Aldrich) and Banana agar medium (BAM)[52] (Supplementary Table 2).

For the liquid culture, BACTEC ™ vials emptied of their contents and filled with Sabouraud broth were used. A total volume of 900 µl of diluted stool sample was inoculated in the liquid culture vials using a syringe needle. Sterile filtered rumen (3 ml) and sheep’s blood (5%) were added to all vials. Then, a total volume of 50 µl of each sample from BD BACTEC vials was cultured on the 5 solid culture media. All culture media were incubated at three different temperatures (22 °C, 28 °C, and 42 °C) and in aerobic/anaerobic conditions.

Fungal identification by MALDI-TOF mass spectrometry

Fungi were identified using Matrix Assisted Laser Desorption/ionization Time-Of-Flight (MALDI-TOF) mass spectrometry (BrukerDaltonics, Germany). A standard method was used for protein extraction, according to the manufacturer’s recommendations[27,29]. One μl of each protein extraction supernatant and the matrix HCCA (a-cyano-4-hydroxycinnamic acid) were spotted on the MSP96 target (Bruker Daltonics, Germany) with a positive control (Candida albicans) and a negative control (the matrix). The data acquisition was performed on the MALDI Biotyper (Microflex LT system; Bruker Daltonics GmbH, Bremen, Germany) using Flex ControlTM software and MALDI BioTyper RTC identification software (Bruker Daltonics).

Molecular fungal identification

To confirm the results of the MALDI-TOF mass spectrometry, direct internal transcribed spacer (ITS) analysis was performed for all fungal species isolates. The ITS1F (CTTGGTCATTTAGAGGAAGTAA) and ITS4R (TCCTCCGCTTATTGATATGC) primers were used for PCR and sequencing. Total DNA was extracted from all strains using a modified protocol EZ1 Advanced DNA bacteria extraction (QIAamp, QiagenInc, Germany). Specifically, 210 µl buffer G2, 10 µl of K proteinase and 100 µl of glass beads (Glass beads, acid-washed, SIGMA ®) were added to an Eppendorf tube of 1.5 ml containing 2 or 3 fungal colonies. An additional heating step of 5 minutes’ duration at 95 °C was performed for filamentous fungi. The supernatant was recovered and we followed the automatic Protocol EZ1 DNA bacteria (QiagenInc). The DNA fragments were amplified with a thermocycler (ThermoCycle2720 Applied Biosystems ®). PCR mix without DNA was used as a negative control for each PCR run. The PCR product (50 µl final volume) contained 5 µl DNA and 45 µl of mix. The mix contained 5 µl of DNA polymerase buffer (Qiagen), 2 µl of MgCL2 (25 mM), 0.25 µl FASTATaq DNA polymerase (1.25 U) (Qiagen), and 1 µl of each primer (Eurogentec, Seraing, Belgium). The PCR program used was as follows: 1 cycle of enzyme activation at 96 °C for 5 minutes, followed by 40 cycles at 96 °C for 30 seconds, 50 °C for 30 s, 72 °C for 1 min and a final cycle at 72 °C for 5 min. The fragment amplification was visualized by electrophoresis on a 1.5% agarose gel. The PCR products were then purified using the Nucleo-Fast® 96 PCR Kit (Marcherey-Nagel, Hoerdt, France) and the purified PCR products were sequenced using the Big Dye® Terminator V1.1 Cycle Sequencing Kit (Applied Biosystems) on a ABI PRISM 3130 automated sequencer (Applied Biosystems®). The sequences were corrected with CodonCode software (6.0.2), and then the sequences were compared with a blast search against the UNITE fungal ITS database (https://unite.ut.ee/).

ITS1 and ITS2 amplicon-based metagenomics

Amplification processes

Total DNAs were extracted using a modification of the Qiagen stool procedure and the Qiamp® DNA Stool Mini Kit (Qiagen, Courtaboeuf, France) as described previously[9]. Both the ITS region 1 and ITS region 2 were amplified with sets of primers (as detailed in Supplementary Table 7) using two different annealing temperatures (52 °C and 55 °C), and three replicated PCRs for each annealing temperature were performed. Genomic amplifications and purifications of PCR products were performed as described previously[9]. All amplicons obtained from each replicated PCR and at both annealing temperature were pooled together for each ITS1 and ITS2 PCR. The total amplicons from both ITS PCRs (ITS region 1 and ITS region 2) were then pooled and sequenced by MiSeq technology.

Illumina MiSeq sequencing

PCR amplified templates from genomic DNA using the ITS primers with overhang adapters (Supplementary Table 7). Samples were amplified for both ITS regions by Taq Phusion (Thermo Fisher Scientific Inc, Waltham, MA U.S.A.) and visualized on the Caliper LabchipII device by a DNA 1 K LabChip. After purification on AMPure beads, the concentrations were measured using high sensitivity Qubit technology (Beckman Coulter Inc, Fullerton, CA, USA), then diluted to 0.2ng/µl. Illumina sequencing adapters, and dual-index barcodes were added to each amplicon using a subsequent limited cycle PCR on 1 ng of each PCR product. The purified libraries were then normalized according to the Nextera XT protocol (Illumina Inc, San Diego, CA, USA). The multiplexed samples were pooled into a single library for sequencing on the MiSeq (Illumina Inc, San Diego, CA, USA). Automated cluster generation and paired-end sequencing with dual index reads was performed in a single 39-hour run in 2 × 250 bp. Total information of 5.7 Gb was obtained from a 588 K/mm2 cluster density with cluster passing quality control filters of 88.5% (11,114,000 clusters). Within this run, the index representations for all the samples were determined between 0.28 and 3.05% with an average of 1.1%. Raw data were configured in fastaq files for R1 and R2 reads between 31 085 to 339 287 paired end reads with an average of 123 741 paired end reads.

Bioinformatics and Statistical analyses

The Bioinformatics analyses were run with minor modifications as described by Balint et al.[53]. PandaSeq was used to assemble the generated paired-end reads by Illumina MiSeq sequencer after trimming both ends[54]. The reads were then reoriented into 5′-3′ directions. Demultiplexing was performed with fqgrep (https://github.com/indraniel/fqgrep). Initial denoising was performed with a 99% similarity clustering with the heuristic clustering algorithm uclust 2.1, implemented in usearch v.6.0.203[55]. De novo chimera detection was performed with the uchime algorithm[56]. OTU picking was performed at 97% sequence similarity with uclust 2.1. ITSx 1.0.11[57] was used to extract fungal ITS1 and ITS2 reads from ITS sequences. The ITS1 and ITS2 sequences were compared using blast search against the UNITE fungal ITS database (https://unite.ut.ee/). We parsed the blast outputs in MEGAN 5.10.5[58] for initial taxonomic screening (minimum reads: 1, minimum score: 170, and upper percentage: 5) and retained OTUs with supported fungal origin for downstream analysis. Alpha-diversity indices, such as Chao-1 and Shannon were calculated by using the PAST3 software package[59]. Statistical analyses were performed using GraphPad Prism version 7.00 for Windows (GraphPad Software).

Nucleotide Sequence Accession Numbers

The read sequences obtained from llumina MiSeq were submitted to the NCBI Sequence Read Archive (SRA) under accession numbers SRR5141534-SRR5141561.