|Home | About | Journals | Submit | Contact Us | Français|
The human intestine is home to a diverse range of bacterial and fungal species, forming an ecological community that contributes to normal physiology and disease susceptibility. Here, the fungal microbiota (mycobiome) in obese and non-obese subjects was characterized using Internal Transcribed Spacer (ITS)-based sequencing. The results demonstrate that obese patients could be discriminated by their specific fungal composition, which also distinguished metabolically “healthy” from “unhealthy” obesity. Clusters according to genus abundance co-segregated with body fatness, fasting triglycerides and HDL-cholesterol. A preliminary link to metabolites such as hexadecanedioic acid, caproic acid and N-acetyl-L-glutamic acid was also found. Mucor racemosus and M. fuscus were the species more represented in non-obese subjects compared to obese counterparts. Interestingly, the decreased relative abundance of the Mucor genus in obese subjects was reversible upon weight loss. Collectively, these findings suggest that manipulation of gut mycobiome communities might be a novel target in the treatment of obesity.
The human gut is a robust ecosystem composed of a dynamic microbial community, collectively termed microbiota, which have important roles in the acquisition of energy from foods and the regulation of host physiology through immune modulation1,2. The finely tuned equilibrium between the host and microbiota may be disrupted when changes in the latter, described as “dysbiosis”, occur in the context of prevalent disorders such as obesity, diabetes and metabolic inflammation2,3.
Bacteria are the most abundant components of the human gut microbiome1,4,5. The composition of distal gut microbiota is known to be altered in obesity and obese individuals generally harbour a decreased ratio of Bacteroidetes to Firmicutes1,2,4. The relationship between gut microbiota and the emergence of obesity4,6,7 might be causal, opening a myriad of possibilities for novel treatment approaches.
The mycobiome, referring principally to the fungal component of microbiota, comprises approximately 0.03–2% of total gut microorganisms and is an integral part of the gastrointestinal tract8,9,10,11. The diversity of fungi inhabiting the human gut remains poorly explored9,12,13,14. A fungal cell is >100-fold larger than a bacterial cell, and thus fungi represent substantially greater biomass than is suggested by the number of available genomes3.
Next generation sequencing has been valuable for characterizing the human gut mycobiome3,12,15,16. Studies have shown that the human gut is home to more than 66 genera and 184 species of fungi, with Candida, Saccharomyces and Cladosporium being particularly common3,9,11,12,15,17,18. In mice, the majority of gut fungi are indigenous to the intestine since only a few of the most common gut fungi were found in mouse food16.
Mycobiome dysbiosis is relevant in inflammatory diseases such Crohn’s disease or ulcerative colitis16,17,18. As yet, no studies have addressed the role of gut fungal microflora in obesity and related diseases.
Here, we used Multitag Pyrosequencing (MTPS) to evaluate fungal diversity in faecal samples from obese and non-obese patients using internal transcribed spacer (ITS) primers, which have broad fungal specificity. We found that the faecal mycobiome is disturbed in obese patients in comparison with non-obese subjects, in parallel with alterations in glucose and lipid metabolism. We observed an optimal grouping of individuals according to genus abundance, that clustered together with metabolic phenotypes. Additionally, we found a dynamic relationship between adiposity and the gut mycobiome after weight loss, indicating that manipulation of gut mycobime communities could be a novel approach in the treatment of obesity.
A total of 52 Caucasian subjects were recruited at the Hospital Universitari Dr. Josep Trueta, Girona (Spain) (Supplementary Table 1). The patients were matched by age and sex and stratified according to BMI in order to evaluate the associations between faecal fungi and obesity. Clinical and analytical characteristics according to obesity status are shown in Table 1.
The Internal Transcribed Spacer (ITS)-based sequencing effort yielded 126704 sequence reads, of which 102891 corresponded to known fungi belonging to 5 different phyla, 13 classes, 50 families and 75 different genera (Supplementary Table 2).
Fungal sequences were detected in every sample analyzed. Less than 1% of sequences corresponding to unknown fungi were detected in 48% of subjects and between 1 and 20% in 17.3% of subjects. The majority of classes, families and genera belonged to the phylum Ascomycota, which corresponds to the largest fungi phylum, and were detected in all samples studied. The phyla Basidiomycota and Zygomycota were detected in 90.38% and 42.30% of samples, respectively. The phyla Chytridiomycota and Neocallimastigomycota, which include a small number of fungi, were present only in two subjects.
At the class level, Saccharomycetes was the most predominant, and was detected in every sample analyzed, followed by Eurotiomycetes and Dothideomycetes, present in 92.3% and 61.53% of the samples, respectively; all belonging to the phylum Ascomycota. The most abundant class of the Basidiomycota phylum was Agaricomycetes, which was present in 46.15% of the samples.
The most abundantly detected families were: Aspergillaceae (86%), belonging to the class Eurotiomycetes, and Dipodascaceae (52%) and Saccharomycetaceae (67.3%) belonging to the class Sacharomycetes. Finally, the family Mucoraceae, one of the three family sequences from the phylum Zygomycota, was identified in 38% of samples.
The most prevalent genera in our analysed samples were Penicillium (present in 73% of the samples), followed by Candida (55%), Saccharomyces (55%), Mucor (38%) and Aspergillus (35%). The relative proportions of the most abundant genera detected by pyrosequencing in each patient are shown in Supplementary Figure 1.
No differences were detected in the richness of the mycobiome between non-obese and obese subjects (Fig. 1a). However, family biodiversity was significantly lower in obese subjects compared with non-obese individuals, and a tendency towards increased biodiversity at other levels was also found in non-obese individuals (Fig. 1b).
Obese and non-obese patients could be distinguished by their specific fungal composition (Fig. 1c,d). The relative abundance of the two major phyla, Ascomycota and Basidiomycota, was not significantly different between obese and non-obese subjects; however, the minor phylum Zygomycota was significantly under-represented among obese subjects (p=0.031) (Fig. 1e).
At the class level, Tremellomycetes appeared only in a subset of obese subjects, whereas it was completely absent in non-obese subjects (p=0.028). In contrast, the class Agaricomycetes was significantly more abundant in non-obese patients compared to obese patients (p=0.017) (Fig. 1f).
At the family level, no significant differences were observed between obese and non-obese subjects. Nonetheless, Aspergillaceae (15.91%) and Mucoraceae (7.58%) were the most prevalent families in non-obese subjects, while Dipodascaceae (7.31%), Aspergillaceae (6.93%) and Saccharomycetaceae (6.13%) were the most abundant in obese subjects (Fig. 1g).
Candida, Nakaseomyces and Penicillium (present in 6.57%, 4.33% and 3.12%, respectively) were the most abundant genera detected in obese patients (Fig. 1h)., Mucor was the most prevalent genus in non-obese patients (8.07%), followed by Candida and Penicillium (5.79% and 2.49% respectively).
We used multivariate analyses to further evaluate the potential differences in the composition of the mycobiome with respect to obesity. A Partial Least Square Discriminant Analysis (PLSDA) revealed that a model with a high accuracy (three component model accuracy 0.77; R2 0.83 and Q2 0.027) could distinguish obese from non-obese mycobiomes (Fig. 1c,d). In order to examine potential overfitting, we also tested as alternative systems for classification random forest-type classifications and recursive support vector machine, by using all genera of the samples. Regarding random forest-type classification, the system produced a relatively low classification error for obese individuals (18%), with high abundance of Aspergillus genus and low abundance of Mucor, Penicillium, Saccharomyces and Eupenicillium genera contributing to obesity (Supplementary Fig. 2). Concerning the recursive support vector machine, this system generated a 10-genus model with an overall error rate near 30% (Supplementary Fig. 3). In this model, low abundance genera, such as Moniliella, Paralepista, Fuscosporia, Limonomyces were present in most of the models tested, but also low Eupenicilium abundance was present in obese individuals. Further, by including a model of only Eupenicillium, Mucor, Penicillium, Moniliella and Eurotium, PLSDA model accuracy was 79.3%, with low discovery rate (permutation test, p=0.04), suggesting that all these genera are fairly associated with obesity (data not shown).
Next, the relative genus abundances were used for optimal grouping (clustering) of individuals according to the Calinksi-Harabasz index. Interestingly, these tests uncovered the presence of 3 potential clusters (Fig. 2a,b) that differed according to the relative abundance of Penicillium spp., Mucor spp., Aspergillaceae, Saccharomycetes, Eurotiomycetes and Zygomycota (Fig. 2c). Strikingly, these clusters were associated with body fatness (BMI, DEXA-measured fat mass and android fat mass), and also with fasting triglycerides and HDL- cholesterol (Fig. 2d). The relationships among anthropometric and metabolic parameters were more marked between subjects harbouring cluster 1 (Fig. 3a).
Heat map analysis demonstrated significant associations among anthropometric and metabolic parameters and fungal taxa (Fig. 3b). The relative abundance of the phylum Zygomycota, class Agaricomycetes, families Mucoraceae and Nectriaceae and genus Mucor and Penicillium correlated negatively with parameters of body fatness such as BMI, fat mass, android fat mass and hip circumference. Conversely, the relative abundance of the classes Sacharomycetes, Tremellomycetes, Cystobasidiomycetes, and families Erythrobasidiaceae and Dipodascaceae and genus Aspergillus was positively associated with adiposity. Moreover, serum total cholesterol, LDL-cholesterol and fasting triglycerides were negatively linked to the relative abundance of sequences belonging to phylum Zygomycota, and also to the class Eurotiomycetes, families Mucoraceae, Hypocraceae and genus Mucor. Conversely, the relative abundance of those sequences belonging to the family Dipodascaceae was positively associated with serum total cholesterol and fasting triglycerides. The relative abundance of fungi belonging to the class Eurotiomycetes, family Aspergillaceae and genus Penicillium was positively correlated with HDL-cholesterol, whereas the abundance of the classes Saccharomycetes, Tremellomycetes, Cystobasidiomycetes and family Erythrobasidiaceae was negatively correlated with this parameter (Fig. 3b).
Fascinatingly, the relative abundance of some fungi was linked to parameters of glucose metabolism. Specifically, the relative abundance of fungi of the class Agaricomycetes, and families Mucoraceae, Ceratocystidaceae, Corticiaceae, Debariomycetaceae, genus Mucor, Monilliela, Ceratocystis and Eupenicillium were negatively associated with HOMA value, glycated haemoglobin and AUC of glucose and insulin during oral glucose tolerance testing. Conversely, the relative abundance of the genus Eurotium was positively associated with fasting glucose and glycated haemoglobin, while the phylum Ascomycota was positively associated with fasting insulin, and the genus Aspergillus correlated positively with AUC insulin. Of note, we also observed associations with parameters of metabolic inflammation (C-reactive protein and lipopolysaccharide binding protein) (Fig. 3b).
Amongst obese individuals, we observed significant associations between the relative abundance of sequences belonging to the Eurotiomycetes class (belonging to Ascomycota phylum) and metabolic parameters (Supplementary Table 3). Interestingly, although the mean relative abundance of Eurotiomycetes was similar between obese (13.02%) and non-obese (13.40%) subjects, 40% of obese patients had an abundance lower than 1%. We observed that obese subjects with less than 1% of sequences belonging to Eurotiomycetes in faecal samples had a more pronounced dyslipidemic profile, increased fasting triglycerides and increased total cholesterol concomitant with increased fasting hyperinsulinemia, compared with obese patients with Eurotiomycetes >1% (p=0.022, p=0.036 and p=0.007, respectively) (Supplementary Table 3) (Fig. 4a,b), and a tendency towards increased HOMA and LDL-cholesterol (p=0.063 and p=0.055, respectively) (Supplementary Table 3). Interestingly, no significant differences in fasting insulin and fasting triglycerides were found between obese subjects with Eurotiomycetes present in >1% and non-obese subjects (Fig. 4a,b).
We performed a plasma metabolomics profile in 8 subjects (Supplementary Table 4). Several metabolites differed significantly in concentrations between individuals with Eurotiomycetes abundance <1% and those with >1% (Fig. 4c). Remarkably, the relative abundance of Eurotiomycetes was associated with plasma concentrations of N-acetyl-L-glutamic acid, caproic acid and hexadecanedioic acid, among other metabolites (Fig. 4d–f). The abundance of Mucor genus was also found significantly correlated with certain metabolites (Fig. 4g).
Since the Mucor genus was significantly more abundant in non-obese subjects, we investigated the potential relationship between gut mycobiome ecology and body fat. We evaluated the relative abundance of Mucor spp. in 14 patients by real time PCR before and after diet-induced weight loss. Results showed that the relative abundance of Mucor spp. for each subject paralleled the degree of weight loss (Fig. 5a). Additionally, pyrosequencing of Mucor positive patients showed that M. fuscus, M. circineloides, M. velotinosus and M. racemosus were the species identified in Mucor genus positive faeces. In non-obese patients, M. racemosus and M. fuscus were significantly more present (p=0.0099 and p=0.0316, respectively) when compared with obese patients (Fig. 5b).
The relative abundance of M. racemosus was negatively associated with BMI, waist circumference, hip circumference, % total fat, fat android distribution, fat gynoid distribution, LDL-cholesterol, fasting triglycerides, uric acid and CRP (Fig. 5c).
This study shows that the mycobiome of obese subjects has an increased presence of the phylum Ascomycota, class Sacharomycetes and families Dipodascaceae and Saccharomycetaceae and, an increased relative abundance of fungi belonging to class Tremellomycetes, compared with non-obese subjects. Although fungi are known to be linked to a number of gastrointestinal diseases11,14,15,17,18, the association between fungal microflora and obesity is novel. A clear tendency towards decreased diversity was observed in obese subjects. These findings are consistent with previously reported data on bacterial diversity in the obese state19 and are contrary to the observed increase in fungal diversity in patients with inflammatory bowel disease or chronic hepatitis B13,17,18. Overall, these findings indicate that the fungal gut mycobiome seems to be disturbed in obese patients, in association with alterations in lipid and glucose metabolism.
The relative abundance of some fungi was linked to adiposity and related metabolic disorders including insulin resistance, dyslipidaemia, blood pressure and inflammatory activity. For example, the phylum Ascomycota, classes Sacharomycetes, Tremellomycetes and Cystobasidiomycetes, families Erythrobasidiaceae and Dipodascaceae and genera Aspergillus, Eurotium and Rhodotorula, all increased with the occurrence of metabolic abnormalities. Conversely, the relative abundance of fungi belonging to the phylum Zygomycota, classes Agaricomycetes and Eurotiomycetes, families Mucoraceae, Nectriaceae, Ceratocystidaceae, Corticiaceae, Debariomycetaceae and Hypocraceae, and genera Mucor, Penicillium, Monilliela and Ceratocystis were associated with protection from these metabolic disorders. The finding that systemic inflammation is associated with obesity suggests that the obesity-related bacterial gut composition has a proinflammatory effect19. This is similar to our observations with Tremellomycetes. This class was not only significantly more abundant in obese subjects, but also correlated positively with inflammatory parameters. Of note, some metabolites from the Penicillium genus have been shown to exhibit anti-inflammatory and insulin sensitizing activities20.
Interestingly, the presence of some fungal communities was associated with metabolically healthy obese subjects. Although the relative abundance of Eurotiomycetes was similar in obese and non-obese subjects, 40% of obese subjects had an abundance less than 1%, and this was associated with a worse metabolic glucose and lipid profile. Moreover, those obese subjects with >1% Eurotiomycetes had similar fasting insulin and fasting triglycerides compared with non-obese subjects.
N-acetyl-L-glutamic acid is known to be metabolized by some members of the Ascomycota phylum, to which Eurotiomycetes belongs21. Interestingly, N-acetyl-L-glutamic acid salts have been shown to exhibit natriuretic effects, resulting in the lowering of blood pressure22. Furthermore, hexadecanedioic acid, associated with a decreased abundance of Eurotiomycetes, has been demonstrated to exhibit antimycotic activity23. Further studies should be performed to experimentally address whether the metabolites described here correlate only with the expansion (or decrease) of specific fungi, or whether they have a direct effect on these fungi.
Candida spp. have been successfully identified from the intestine of healthy individuals9,17. Some studies have suggested a link between expansion in Candida spp. and diabetes24 and inflammatory disorders of the gastrointestinal tract17,25, with a possible active role for Candida albicans26. In this sense, we have not observed differences of the relative expression of Candida albicans in obesity (data not shown).
In our study, the Mucor genus was the most prevalent genus in non-obese subjects; specifically, M. racemosus and M. fuscus were the more representative species in these subjects. Mucor sp. has a cell wall polysaccharide composition based on chitin-chitosan27, and multi-functional roles of this polysaccharide as a protective agent against obesity have been described28. M racemosus has been described as source of chitosan29. Future animal studies will help to evaluate the role of the identified Mucor spp. in obesity.
Mechanistically, we also demonstrate that the relative abundance of the Mucor genus increased after weight loss in obese subjects in a manner analogous to Bacteroidetes, which increases also after weight loss4,30. However, a larger study should be carried out to evaluate the significance of these findings.
At the bacterial level, three predominant enterotypes have been identified based on variations in the levels of three different genera (Bacteroides, Prevotella and Ruminococcus)10,31. Here we identified 3 different clusters; strikingly, cluster 1 was significantly associated with adiposity markers and dyslipidemia, hinting at a possible target for treatment. We observed that the major shift in fungal populations associated with cluster 3 mirror similar dramatic changes in the bacterial communities (ratio Firmicutes/Bateriodetes) than cluster 1 and 2 (data not shown). One could argue that 4 samples that fall in the less abundant cluster 3 are actually outliers. In order to exclude that proposed clusters were not derived from outlying measures, we performed Grubb’s outlier test on the matrix of Jensen-Shannon distances, by employing the R package “outliers”, and found no significant outliers. However the number of subjects studied within this cluster was too low to reach to any conclusion. We acknowledge that further cluster validation is required for a generalization of these results as potential “enteromycotypes”.
Finally, fungal-fungal interaction analysis showed that, in obese patients, the Ascomycota phylum negatively correlated with Basidiomycota and Zygomycota phyla. Noteworthy, among non-obese subjects, a significant negative association was observed between Pichiaceae and Dipodascaceae. The same antagonistic relationship between Candida (Family Dipodascaceae) and Pichia (Family Pichiaceae) has been observed in oral mycobiota of HIV-infected patients32.
In summary, this study represents the first analysis of the human mycobiome in obesity and associated metabolic disorders. There is limited knowledge on gut bacterial-fungal interactions and their role in health and disease. An understanding of cross-kingdom microbial interactions within the context of health and disease holds considerable promise to facilitate the discovery of potential preventative and therapeutic targets. Further analysis on the metabolic road-maps involving fungal and bacterial genomes will be needed to unequivocally unravel new avenues in inter-kingdom metabolic dependencies.
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.4kg/m2. Subjects were instructed to maintain a 20kcal/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.3kg/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.
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 5minutes of rest. Three readings were obtained and the mean value was used in the analyses.
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 ≥30kg/m2. Android fat mass was calculated using regions of interest at the abdominal level in DEXA scans.
After 8h 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 10mU/l and 3.4% at 130mU/l. The inter-assay CVs were 6.9 and 4.5% at 14 and 89mU/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 75g oral glucose tolerance test was performed after an overnight fast and venous blood samples were drawn at time points 0 30, 60, 90 and 120min 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 75g 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.02mg/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 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 2ml 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 40seconds 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).
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 GCTRYGTTCTTCATCGDTGC. Primer set B; H2SeqFb GCA TCG ATG AAG AACRYA 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 15min followed by 35 cycles of 95° for 20sec, 62° for 30sec and 72° for 60sec, and a final step of 72° for 10min. 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.
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 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.
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: 2min at 50°C for Uracil-DNA Glycosylase treatment, 10min at 95°C for Taq activation, 15s at 95°C for denaturation and 1min 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). Ct 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).
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×300bp. 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.
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 1h at −20°C and centrifuged for 3min at 12000g. 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 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.
How to cite this article: Mar Rodríguez, M. et al. Obesity changes the human gut mycobiome. Sci. Rep. 5, 14600; doi: 10.1038/srep14600 (2015).
We wish to acknowledge the kind collaboration of the people who voluntarily participated in the study. This study was supported by projects from the “Fondo de Investigación Sanitaria” (FIS): PI 11/00701; PI 11/00049; PI 14/00465 (M.R.CH), co-financed by the European Regional Development Fund (ERDF). MRCH is supported by the Research Stabilization Programme of the Instituto de Salud Carlos III (ISCIII) co-financed by Institut Català de Salut (ICS) in Catalonia. M del Mar Rodriguez is supported by a “Juan de la Cierva” Fellowship JCI-2011-11488.
Author Contributions The experiments were conceived and designed by J.M.F.R., M.R.C.H. and F.J.C. The experiments were performed by: M.M.R. and M.J. The data was analyzed by: D.P., P.M.G., F.J.C., M.M.R. and M.P.O. The following authors contributed in buying reagents or materials or provided analysis tools: M.R.C.H., J.M.F.R., F.J.C., J.V., M.J., R.P., G.X., W.R. and E.E. Finally, M.M.R., M.R.C.H. and J.M.F.R. where responsible for writing the paper.