Subjects and Sample Collection

Two cohorts were used to perform the study. In cohort 1, 52 subjects were recruited at Hospital Universitari Dr. Josep Trueta, Girona (Spain) (Supplementary Table 1). All subjects were Caucasian and reported that their body weight had been stable for at least 3 months before the study. These subjects had no systemic disease other than obesity or dyslipidaemia, were free of any infection in the month before the study and did not undergo treatment with drugs that affect glucose metabolism or antibiotics. Liver and renal diseases were specifically excluded by biochemical work-up.

Cohort 2 included 14 overweight and obese patients from the Hospital Universitari Dr. Josep Trueta. The subjects were 7 men and 7 women, aged 54.2 ± 8.7 years, with a mean BMI of 33.4 ± 7.4 kg/m2. Subjects were instructed to maintain a 20 kcal/kg balanced diet. Measurements and samples were taken at baseline and 16 weeks after starting the diet. At the end of this period, mean BMI was 31.5 ± 7.3 kg/m2 (p = 0.005 vs. baseline).

All experiments were performed in accordance with approved guidelines and regulations. All experimental protocols were approved by the Ethics Committee of the Hospital Universitari Dr. Josep Trueta (Comitè d’Ètica d’Investigació Clínica, CEIC, ethical approval number 2009046). Informed consent was obtained from all participants.

Clinical and biochemical variables

Each patient underwent evaluation of anthropometric and laboratory parameters on the same day. Blood pressure was measured using a blood pressure monitor (Hem-703C, Omron, Barcelona, Spain), with the subject seated and after 5 minutes of rest. Three readings were obtained and the mean value was used in the analyses.

Body composition

Fat mass was evaluated by dual-energy x-ray absorptiometry (DEXA) using a GE Lunar Prodigy Oracle densitometer (GE Healthcare, enCore software version 13.2). Whole body composition (fat mass and fat-free soft tissue mass) was obtained by trained personnel according to standard procedures. Obesity was considered as a BMI ≥30 kg/m2. Android fat mass was calculated using regions of interest at the abdominal level in DEXA scans.

Analytical methods

After 8 h of fasting, blood was drawn for the measurement of fasting plasma lipids, glucose and insulin. Serum glucose concentrations were measured in duplicate by the glucose oxidase method using a Beckman Glucose Analyser II (Beckman Instruments, Brea, CA). Intra-assay and inter-assay coefficients of variation (CV) were less than 4%. Plasma lipid determinations were performed on a Roche/Hitachi Cobas c 711 autoanalyzer (Roche Diagnostics GmbH, Mannheim, Germany). Total serum cholesterol was measured by an enzymatic colorimetric cholesterol esterase/cholesterol oxidase/peroxidase reaction (Cobas CHOL2). HDL cholesterol was quantified by a homogeneous enzymatic colorimetric esterase/cholesterol oxidase/peroxidase reaction (Cobas HDLC3). LDL cholesterol was calculated using the Friedewald formula. Total serum triglycerides were measured by an enzymatic colorimetric glycerol phosphate oxidase and peroxidase reaction (Cobas TRIGL).

Glycated haemoglobin (HbA1c) was measured by high-pressure liquid chromatography using a fully automated glycosylated haemoglobin analyzer system (Hitachi L-9100). Serum insulin was measured in duplicate using a monoclonal immunoradiometric assay (Medgenix Diagnostics, Fleunes, Belgium). The intra-assay CV was 5.2% at a concentration of 10 mU/l and 3.4% at 130 mU/l. The inter-assay CVs were 6.9 and 4.5% at 14 and 89 mU/l, respectively. The homeostatic model assessment insulin resistance index (HOMA-IR) was calculated by the formula: (serum glucose (mmol/l) x serum insulin (mU/l))/22.5. A standard 75 g oral glucose tolerance test was performed after an overnight fast and venous blood samples were drawn at time points 0 30, 60, 90 and 120 min for determination of plasma glucose and insulin. Area under the curve of glucose (AUC-glucose) and insulin (AUC-insulin) concentrations during the 120 min of the 75 g oral glucose tolerance test was then calculated using the trapezoid method. C-reactive protein (ultrasensitive assay; Beckman, Fullerton, CA) was determined by a routine laboratory test, with intra- and interassay CVs of 4%. The lower limit of detection is 0.02 mg/l. Serum ferritin was measured by microparticle enzyme immunoassay (AxSYM™; Abbot Laboratories) with intra- and interassay CVs < 6%. Serum lipopolysaccharide binding protein (LBP) levels were measured with an ELISA kit (HK315-02, HyCult Biotech Inc., Uden, The Netherlands). Serum samples were diluted and assayed according to the manufacturer’s instructions. Intra- and inter-assay CVs for all the determinations were between 5 and 10%. Serum alanine aminotransferase (ALT), aspartate aminotransferase (AST) and gamma-glutamyltransferase (GGT) levels were determined using enzymatic methods. Uric acid was determined by routine laboratory tests.

Stool processing and DNA extraction

Stool specimens were obtained from patients using sterile containers and were immediately frozen in liquid nitrogen and stored at −80 °C until analysis. Samples were processed individually using the Fast DNA Spin Kit for faeces (MP Biomedicals, Solon, OH). Briefly, a frozen aliquot (400mg) of each sample was added to a 2 ml tube containing 825 μl sodium phosphate buffer, 275 μl of pre-lysis solution and Lysing matrix E, a mixture of ceramic and silica particles designed to efficiently lyse all stool microorganisms. Each extraction tube was agitated twice for 40 seconds using a Fast Prep FP120 instrument at a speed setting of 6 to ensure proper extraction of fungal DNA, a crucial point in the methodology. Tubes were cooled on ice between the different agitation procedures. DNA extraction was then carried out following the manufacturer’s instructions. The quantity and quality of isolated DNA was determined with a Nanodrop ND-1000 spectrophotometer (Nanodrop Technologies, Wilmington, DE, USA).

Pyrosequencing analysis

Fungal genomic sequences present in faecal samples were identified by amplification of the internal transcribed spacer-based region (ITS). The ITS primers used span the region between the 3′ end of the 18S gene, including the entire 5.8S gene and end of the 5′ region of 28S gene. PCR reactions were performed in 15 μl reaction volumes containing 0.33 μM of primer. Two different PCR reactions were performed per sample, using two different pair of primers in a 384-well plate. Both sets of primers contained a universal tag sequence in the second PCR step (5′ = AGGTCAGGATCAACGCTCAAG, 3′ = CATCTTGCATGATCCAACCTTC). Primer set A; H1SeqF GTCATTTAGAGGAAGTAAAAGTCGTAACAAGG and H1SeqRb GCT RY GTTCTTCATCG D TGC. Primer set B; H2SeqFb GCA TCG ATG AAG AAC RY A GC and primerH2SeqRb TTC TTT TCC TCC GCT TAT TGA TAT GC. Standard PCR cycling was performed on a Veriti 384-well Thermal Cycler (Applied Biosystems) with an initial step at 95° for 15 min followed by 35 cycles of 95° for 20 sec, 62° for 30 sec and 72° for 60 sec and a final step of 72° for 10 min. PCR products were visualized using a Qiaxcel system (Qiagen).

PCR products from the first step were diluted 1/5 with water and were used as templates in the second PCR step in order to tag every PCR product with a sample specific tag multiplex identifier (MID) consisting of a 10 base pair specific sequence. Two μl of the diluted PCR product were used for the second PCR reaction, containing 0.3 μM of each MID. Standard PCR cycling was performed as before. PCR products were visualized using a Qiaxcel system, purified using size-exclusion AMPure SPRI DNA-binding paramagnetic beads (Agencourt Bioscience Corp. Beckman Coulter) and quantified in 96 well-format with the QuantiFluor-ST Fluorometer (Promega) using a PicoGreen® assay (InvitroGen). Samples were then pooled at an approximate equimolar concentration. After pooling, the library was further purified with the Pippin Prep size-exclusion system (SageScience), quantified in 96-well format with the QuantiFluor-ST Fluorometer (Promega) and diluted to 5 × 105 molecules/μl. Five μl of the pools were then subjected to emulsion PCR (emPCR) using the 454 Junior Titanium Series LibA emPCR kit (Roche diagnostics) to perform last amplification step.

Data Analysis

The 454 multitag pyrosequencing (MTPS) data was de-multiplexed by sorting the sequences into bins based on the barcodes in the samples. A quality control was carried out on the reads for each sample using FastQC (v. 0.10.1) to assess potential problems in the data. TagCleaner (v. 0.12) was used to remove the additional tag sequences (sequencing primers and adaptors) that appeared on the raw reads, which could contain deletions or insertions due to sequencing limitations. This step was followed by the trimming of low quality read ends (both 5′ and 3′) with PRINSEQ (v. 0.19.5), removal of short reads, removal of reads with low mean quality and removal of low complexity reads. The cleaned pyrosequencing data was BLASTed against the NCBI nucleotide collection (nr/nt) database using BLAST from the standalone BLAST+ package. We found that 18.8% of the data did not match to fungal sequences. Fungal reads were assigned taxonomically with MEGAN (v. 5.2.3). This step was optimized to maximize the number of reads (Min Support = 1, Min Score = 50.0, Max Expected = 0.01, Top Percent = 5.0, Min Complexity = 0.44). The results were summarized at a variety of taxonomic ranks according to the NCBI taxonomy using the Lowest Common Ancestor (LCA) algorithm in MEGAN. The LCA algorithm assigns taxa to the lowest possible taxonomic rank that presumably reflects the level of sequence variation present in the query sequence compared with reference sequences.

MEGAN classifies each read depending on 4 different types of scenarios: 1) Those reads with a clear species or genera classification get assigned to this taxon. 2) Those with an ambiguous classification (reads match with several different taxa): we report the upper common level of the main matches. 3) Reads in all fungal classes: classified as unknown fungi 4) Those with no matches to fungal sequences: removed from the data analysis but taken into account in the quality control of sequences.

Although parameters were optimized to reach species level as much as possible, when BLAST results classified the sequence as one specific species in most of the hits, MEGAN stops at the higher taxon. We are aware that MEGAN is more conservative than other curated databases, however it is robust enough to be valid for this project and as we are aware of its limitations. Therefore, we decided to provide sequences only up to genus level.

Richness of fungal community was calculated as the number of different taxa for each patient. Biodiversity index was assessed using the Shannon-Weaver index since this is quite insensitive to sample size18. This index was calculated with R-bioconductor package vegan (v. 2.0–10).

Samples clusters according to genus abundance

Samples were clustered based on relative genus abundances using Jensen-Shannon distance and the Partitioning Around Medoids clustering algorithm, previously used for enterotype determination10. The results were assessed for the optimal number of clusters using the Calinski-Harabasz (CH) Index33. These studies, as well as those between class and principal coordinate analyses, were performed by using the r scripts present in http://enterotype.embl.de/enterotypes.html.

Relative quantification of Mucor species by Real-time PCR

Mucor species quantification was performed on a 7900HT Fast Real-Time PCR System using the primers and probes described previously34, with slight modifications. Mucor (F) 5′ GTC TTT GAA CGC AAC TTG CG 3′, Mucor (R) 5′ CCG CCT GAT TTC AGA TCA AAT T 3′ and Mucor probe: 5′ TTCCAATGAGCACGCCTGTT-MGBNFQ 3′. DNA from Mucor circinelloides (CBS277.49), belonging to the mold collection of the Fungal Biodiversity Center (CBS) was used as a positive control. For Mucor detection, PCR was performed in a final volume of 20 μl containing 0.5 μM of each primer of the Mucormycetes species, 0.25 μM of the internal control primers and 0.2 μM of Mucor spp. probe. For 18S (reference gene) amplification, we used primers, probes and conditions described previously35. PCR conditions for Mucor spp., 18S and internal control were: 2 min at 50 °C for Uracil-DNA Glycosylase treatment, 10 min at 95 °C for Taq activation, 15 s at 95 °C for denaturation and 1 min at 61.2 °C for annealing and extension (50 cycles). SDS software 2.3 and RQ Manager 1.2 (Applied Biosystems) were used to analyze the results with the comparative threshold cycle (Ct) method (2−∆∆Ct). C t values for each sample were normalized with the 18S reference gene. All data were expressed as an n-fold difference relative to a calibrator (a mix of DNA from 4 human faecal samples).

Identification of Mucor species by pyrosequencing

For each sample sequenced we created a custom DNA sequencing library specific for Mucor sequences. We performed high-throughput shotgun DNA sequencing using the Illumina MiSeq platform 2 × 300 bp. The data analysis consists of two major stages, the denoising and chimera detection stage and the microbial diversity analysis stage. During the denoising and chimera detection stage, denoising is performed using various techniques to remove short sequences, singleton sequences and noisy reads.

With the erroneous reads removed, chimera detection is performed to aid in the removal of chimeric sequences. Lastly, remaining sequences are then corrected base by base to help remove noise from within each sequence. For this procedure, we used the USEARCH global alignment algorithm and the UCHIME chimera detection software. During the diversity analysis stage, for each sample we determine the taxonomic information for each constituent read. For this procedure we used USEARCH global alignment and UPARSE algorithms and the RDP classifier against a database of high quality sequences derived from the NCBI database. The OTU information was obtained by using MUSCLE aligner and FastTree software.

Metabolomic analysis

Metabolites were extracted from plasma samples of 8 random samples from cohort 1 with methanol according to a previously described method36. Briefly, 30 μl of cold methanol was added to 90 μl of plasma, incubated 1 h at −20 °C and centrifuged for 3 min at 12000 g. The supernatants were recovered, evaporated using a Speed Vac (Thermo Fisher Scientific, Barcelona, Spain) and resuspended in water. We used an ultra-high pressure liquid chromatography (Agilent 1290 LC) system coupled to an electrospray-ionization quadrupole time of flight mass spectrometer (Q-TOF) 6520 instrument (Agilent Technologies, Barcelona, Spain). A column with 1.8 micron particle size was employed and we performed identification of metabolites using the PCDL database from Agilent (Agilent Technologies, Barcelona, Spain), which uses retention times in a standardized chromatographic system as an orthogonal searchable parameter to complement accurate mass data, according to previously published works37.

Statistical analysis

Statistical analysis was performed using version 19 of the Statistical Package for the Social Sciences (SPSS, Chicago, IL). For clinical and anthropometrical variables, data are expressed as mean ± SD. Before statistical analysis, normal distribution and homogeneity of the variances was evaluated using Levene’s test and then variables were log-transformed when necessary. Differences between non-obese and obese subjects were tested with the Mann-Whitney U test for non-normally distributed data. For comparison among clusters and fungi, we used the Mann-Whitney U test, with Bonferroni corrections for multiple comparisons. Differences between cluster 1 and cluster 2 plus 3, in relation to metabolic parameters, were tested using unpaired t-test. Spearman’s correlation coefficient was used to analyze the association between fungi and clinical or metabolic parameters. Partial least square discriminant-analyses were performed using the Metaboanalyst v 3.0 platform37. Two-way ANOVA followed by Bonferroni post hoc comparison test (Prism 5 software; Graph- Pad Software, Inc., San Diego, CA) was used to assess differences in metabolic parameters according to Eurotiomycetes class relative abundance. P < 0.05 was considered significant.