Fecal fungal microbiota alterations associated with clinical phenotypes in Crohn’s disease in southwest China

Study population

Patients in West China Hospital were enrolled from December 2018 to October 2019. Patients with CD were diagnosed by ECCO-ESGAR Guideline for Diagnostic Assessment in CD (Maaser et al., 2019). Those with infectious diseases were excluded, including gastrointestinal or extra-gastrointestinal infections based on clinical symptoms and serological and radiographic results. In addition, patients with cancer, other gastrointestinal disorders, critical illnesses, or intake of antibiotics and/or probiotics in the last 2 months were excluded. Healthy controls (HCs) were volunteers without any suspicious gastrointestinal discomfort, autoimmune diseases, chronic diseases, or history of any kind of medication within 3 months. All included objects were living in southwest China and were in accordance with a spicy dietary habit. One stool sample was collected from each study participant, keeping at ordinary temperature for less than 10 min before being transferred to an ultra-low temperature freezer (–80 °C) until being sent to the Biozeron (China) for DNA extraction and Illumina PE250 sequencing.

Detailed data, including basic information, clinical manifestation, Crohn’s disease activity index (CDAI), Montreal classification (Satsangi et al., 2006), past history of abdominal surgery, and medication, were collected for further analysis. According to the Montreal classification (Satsangi et al., 2006), the biological behaviors of CD were classified into type B1 (non-stricturing, non-penetrating), type B2 (stricturing), and type B3 (penetrating). As types, B2 and B3 often coexist, the biological behaviors of CD patients were divided into two subgroups in this study, including type B1 and type non-B1 (B2 and/or B3). Patients with CDAI score ≥ 150 were enrolled in the group of CD-flare, or else in the group of CD-remission. Analysis was conducted between intestinal fungi from different groups as follows: 1. CD vs. HC (n1 = 45, n2 = 17); 2. CD-flare vs. CD-remission (n1 = 34, n2 = 11); 3. type B1 vs. type non-B1 (n1 = 20, n2 = 25); 4. CD with perianal lesions vs. CD without perianal lesions (n1 = 25, n2 = 20).

Our ethical approval was obtained from the Biomedical Ethics Committee of West China, Sichuan University (NO. 2016-298). Patients were required to give written consents to the study. For full disclosure, the details of the study are published online on the page of the Biomedical Ethics Committee of West China Hospital, Sichuan University.

DNA extraction and Illumina PE250 sequencing

Microbial DNA was extracted using the HiPure Stool DNA Kits (Magen, Guangzhou, China) according to the manufacturer’s protocols. The fecal samples were homogenized in a lysate and further lysed in a high-temperature water bath, releasing DNA into the lysate. The fecal impurities were removed by centrifugation, and the DNA was further purified through the column. Finally, the DNA was eluted by Buffer AE. After monitoring the concentration and purity of extracted DNA on 1% agarose gels, DNA was diluted to 1 ng/μl using sterile water. The fungal microbiota was identified and analyzed by sequencing the internal transcribed spacer (ITS) fragment, which was the target region ITS2 amplified with universal primers: forward primer (5′-GCATCGATGAAGAACGCAGC-3′) and reverse primer (5′-TCCTCCGCTTATTGATATGC-3′) (Bellemain et al., 2010). We performed PCR in a 30 μl solution with 0.2 μm forward and reverse primers and 10 ng template DNA, which contained 15 μL Phusion® High-Fidelity PCR Master Mix (New England Biolabs, Ipswich, MA, USA). The thermocycling included initial denaturation at 98 °C for 1 min, 30 cycles of denaturation at 98 °C for 10 s, primer annealing at 50 °C for 30 s, extension at 72 °C for 5 min, and final extension at 72 °C for 5 min. For qualification and quantification of our PCR products, we mixed PCR products with the same volume of 1× loading buffer containing SYBR green, performed electrophoresis on a 2% agarose gel, and picked samples with the main band brightness between 400 and 450 bp for further experiments. Eventually, we mixed all obtained PCR products in an isopycnic ratio and refined them with the Gene JET gel extraction kit (Thermo Scientific, Waltham, MA, USA).

The sequencing library was created by the NEB Next® Ultra™ DNA Library Prep Kit for Illumina (NEB, USA), complying with the manufacturer’s instructions, adding index codes. The library was estimated on the Qubit@ 2.0 Fluorometer (Thermo Scientific, Waltham, MA, USA) and Agilent Bioanalyzer 2100 system and sequenced on the Illumina MiSeq platform, finally generating 250/300 bp paired-end reads.

Fungal ITS2 rDNA sequence analysis

Paired-end reads were merged as raw tags using FLASH (version 1.2.11) (Magoč & Salzberg, 2011). Paired-end reads were assigned to each sample based on the unique barcodes. Noisy sequences of raw tags were filtered under specific filtering conditions to obtain high-quality clean tags (Bokulich et al., 2013). The software package UPARSE (version 9.2.64) carried out sequence analyses (Edgar, 2013). Conventionally, operational taxonomic units (OTUs) were defined as sequences with ≥97% similarity. After choosing typical sequences for each OTU, we footnoted category data for them by the RDP classifier (version 2.2) (Wang et al., 2007) based on the ITS2 database (version update_2015), with the confidence threshold value of 0.8 (Ankenbrand et al., 2015).

Statistical analysis

Chao1, ACE, Shannon, Simpson, and Good’s coverage index were calculated in QIIME (version 1.9.1) (Caporaso et al., 2010). PD-whole tree index was calculated in Picante (version 1.8.2) (Kembel et al., 2010). All above demonstrate the α diversity between different groups. Principal coordinate analysis (PCoA) based on weighted and unweighted UniFrac distances was generated in the R project Vegan package (version 2.5.3) demonstrating the β diversity between different groups (Oksanen et al., 2010). Adonis (also called Permanova) and Anosim test was calculated in the R project Vegan package (version 2.5.3) (Oksanen et al., 2010). SparCC network analysis visualized the complex cooperative and competitive relationships of fungi within different groups (Friedman & Alm, 2012). Correlation analysis was calculated in the R project psych package (version 1.8.4), ascertaining the association between intestinal fungi and multiple indicators (Revelle, 2022). Differences with a P-value < 0.05 were considered significant. All authors had access to the study data and reviewed and approved the final manuscript.

Study information and patient characteristics

The information of the 45 enrolled patients is described in Table 1 and Table S1. A total of 17 HCs were recruited during the same period, with an average age of 27.6 ± 9.7 years and a male/female ratio of 8/9. There was no significant difference in sex, mean age, Montreal classification, perianal complications, or surgical history between the CD-flare and CD-remission groups.

Altered intestinal fungal microbiota in patients with CD compared to HCs

At the phylum level, there was no significant difference in the composition of Ascomycota, Basidiomycota, and Mucoromycota between the CD and HC groups (P > 0.05). The phyla Rozellomycota and Mortierellomycota were observed in the HC group, while these two phyla were not present in the CD group (Fig. 1A). The relative abundances of Saccharomyces, Clonostachys, and Exophiala in CD patients were higher than those in HCs, while the relative abundance of Candida was approximately half of that in HCs at the genus level (Fig. 1B). At the species level, the relative abundances of Exophiala dermatitidis, Candida tropicalis, Humicola grisea and some unclassified fungal species were higher in CD patients than in HCs, while the relative abundance of Candida albicans was lower than that in HCs (Fig. 1C) (P < 0.05).

Welch’s t-test was used to analyze the mean relative abundance of fecal fungi between these two groups (Fig. 1D). The top 15 genera with significant differences depending on relative abundance were Saccharomyces (P = 0.04), Clonostachys (P < 0.01), Exophiala (P = 0.02), Humicola (P < 0.01), Thermoascus (P < 0.01), Dipodascus (P < 0.01), Cutaneotrichosporon (P < 0.01), Tausonia (P < 0.01), Lophiostoma (P < 0.01), Mucor (P < 0.01), Trichosporon (P < 0.01), Thermomyces (P = 0.04), Cylindrocarpon (P < 0.01), Vanrija (P < 0.01) and Mortierella (P = 0.01). Among them, the top five fungi, depending on relative abundance, all belonged to the phylum Ascomycota. Furthermore, the genera Humicola, Lophiostoma, and Cylindrocarpon only existed in CD patients, all of which were in the phylum Ascomycota. The genus Mortierella was only observed in HC group, with a significant difference compared to CD (P = 0.01). At the species level, there were higher relative abundances of Exophiala dermatitidis (P = 0.02), Humicola grisea (P < 0.01), Thermoascus aurantiacus (P < 0.01), Candida sake, (P < 0.01), Dipodascus australiensis (P < 0.01), Tausonia pullulans (P < 0.01), Cutaneotrichosporon curvatus, (P < 0.01), Mucor racemosus (P < 0.01), Aspergillus sydowii (P < 0.03), and Vanrija fragicola (P = 0.01) in CD patients compared to HCs (Fig. 1E). Among them, the top five fungal species all belonged to the phylum Ascomycota. Humicola grisea, a species belonging to the phylum Ascomycota, was only found in patients with CD. Although two phyla were not present in the CD group, fungi at the genus and species levels belonging to the phylum Ascomycota showed higher relative abundance in the CD group compared with the HC group.

Through the analysis of α diversity indexes, including observed species, Chao1, Shannon, Simpson, and Good’s coverage index, although there was no significant difference between these two groups, patients with CD featured lower fecal fungal richness (observed species, P = 0.27) (Table S2). Principal component analysis (PCoA) based on unweighted UniFrac distance showed the difference in fungal microflora structure between these two groups, indicating that intestinal fungi in some HCs overlapped with those in some CD patients, while some of them were significantly different from the overlapping fungal microflora (Adonis test, P < 0.01) (Fig. 1F). At the genus level, SparCC analysis showed that the number of inter-fungal cooperative or competitive correlations in the CD group was notably reduced compared with that in the HC group (Fig 1G, Table S3). The interacting fungi in the CD group were mainly (7/8) fungi significantly enriched in CD (Figs. 1D and 1G). Among them, the genera Clonostachys, Lophiostoma, and Mucor were positively correlated (Figs. 1D and 1G). The genus Tausonia played a linkage role, was negatively correlated with the above three fungi (absolute value of SparCC correlation > 0.5, P < 0.05), and was positively correlated with the genus Cutaneotrichosporon (SparCC correlation value = 0.51, P < 0.01) (Figs. 1D and 1G).

Fecal fungi in the CD-flare group differed from those in the CD-remission group

The phylum Ascomycota dominated fungi in both CD patients in the active stage (CD-flare, n = 34) and CD patients in the remission stage (CD-remission, n = 11) (Fig. 2A). At the genus level, Candida was proportionable in both the CD-flare group and the CD-remission group. At the same time, Exophiala and Saccharomyces were enriched in the CD-flare group, and Aspergillus and Clonostachys were enriched in the CD-remission group (Fig. 2B). Species Exophiala dermatitidis was also with higher relative abundance in the CD-flare group compared to CD-remission group at the taxonomic levels of its order, family, genus, species, and OTU (P <0.05) (Figs. 2B–2E, Table S4).

Although no significant difference was observed, the observed species richness of the CD-flare was lower than that of the CD-remission (P = 0.08) (Fig. 2F, Table S5). The PD-whole tree index in CD-flare was significantly lower than that in CD-remission (P = 0.03) (Fig. 2G, Table S5). Principal component analysis (PCoA) based on unweighted UniFrac distance showed a significant difference in β diversity between CD-flare and CD-remission (Adonism test, P = 0.027) (Fig. 2H). Furthermore, there was an apparent separation based on weighted UniFrac distance in the fungal community between CD and a part of HC, while another part of HC overlapped with CD, close to the CD-remission fungal community (Adonism test, P = 0.001) (Fig. 2I). All of the above results suggested a progressive change in the fungal community in HC, CD remission, and CD flare.

Correlation analysis between intestinal fungi and laboratory indicators was conducted, including erythrocyte sedimentation rate (ESR), C-reactive protein (CRP), hematocrit (HCT), white blood cell count (WBC), neutrophil count (N), lymphocyte count (N), neutrophil-to-lymphocyte ratio (N/L), platelet count (PLT), lymphocyte percentage, the ratio of PLT to lymphocyte percent, and lymphocyte percent. The results showed that Exophiala dermatitidis was significantly correlated with increased platelet count at the taxonomic levels of its order (Chaetothyriales, R = 0.34, P = 0.02), family (Herpotrichiellace, R = 0.36, P = 0.02), genus (Exophiala, R = 0.38, P = 0.01) and species (Exophiala dermatitidis, R = 0.37, P = 0.01) (Spearman, all P < 0.05) (Fig. 3A). In addition, at the genus level, a positive correlation was observed between Nigrospora (Spearman, R = 0.36, P = 0.02) and ESR. Verticillium was negatively correlated with ESR (Spearman, R = −0.31, P = 0.04) and CRP (Spearman, R = −0.30, P = 0.04) (Fig. 3A). Lophiostoma, which was enriched in the CD group, was positively correlated with the N/L ratio (Spearman, R = 0.32, P = 0.03) (Figs. 1D and 3A). At the species level, Candida albicans, Verticillium dahliae, Wallemia Canadensis, Aspergillus penicillioides, and Nigrospora oryzae correlated with laboratory indicators to some extent (Spearman, P < 0.05) (Fig. 3B). Exophiala dermatitidis was positively correlated with CDAI and laboratory activity index (platelet count) at multiple taxonomic levels, demonstrating that it may be closely related to CD activity.

Differed fecal fungi relating to concerned biological behaviors in CD

In terms of complications, at the genus level, Clonostachys (R = 0.34, P = 0.02), Lophiostoma (R = 0.39, P < 0.01), and Fusarium (R = 0.42, P < 0.01) showed positive correlations with fistula formation, abdominal abscess, and arthritis/arthralgia (Spearman, all P < 0.05) (Fig. 3A). Aspergillus (R = −0.39, P < 0.01), Tausonia (R = −0.39, P = 0.04), Trichosporon (R = −0.33, P = 0.02), Wallemia (R = −0.38, P < 0.01), Holtermanniella (R = −0.30, P = 0.04), and Wickerhamomyces (R = −0.31, P = 0.04) were negatively correlated with the formation of fistula, while Vanrija (R = −0.30, P = 0.047) was negatively correlated with the formation of perianal lesions (Spearman, all P < 0.05) (Fig. 3A). The complex relationships between fungi and complications prompted us to pay more attention to the relationships between fungi and biological behavior and perianal lesions in CD.

Notably, type B1 (n = 20) and type non-B1 (B2/B3, n = 25) were concerned biological behaviors in clinical practice. There was no significant difference in baseline parameters between these two groups (P > 0.05) (Table S6). No significant difference was observed in α or β diversity between these two groups. No remarkable difference was found at the phylum level between these two groups (Fig. 4A). At the genus level, the composition proportions of Saccharomyces and Clonostachys were similar between these two groups (Fig. 4B). At the same time, Candida was significantly enriched in type non-B1 CD compared with type B1 (Welch’s t test, P < 0.05) (Figs. 4B and 4D). At the species level, we did not observe a significant difference in Candida albicans, Candida glabrica, or Candida tropicalis in non-B1-type CD (P > 0.05) (Fig. 4C).

Further analysis was performed grouped by CD patients with perianal lesions (n = 25) and without perianal lesions (n = 20). There was no significant difference in baseline parameters between these two groups (P > 0.05) (Table S7). Fungi at the phylum level were in line with each other (Fig. 4E). The proportion of the genus Candida was similar in these two groups. Clonostachys was enriched in CD patients with perianal lesions, while Saccharomyces and Dipodascus were enriched in CD patients without perianal lesions (Fig. 4F). Furthermore, the fungi significantly enriched in CD patients with perianal lesions at the genus level were Clonostachys (P = 0.02), Humicola (P = 0.02), Lophiostoma (P = 0.03), Fusarium (P = 0.02), Lecanicillium (P < 0.01), and Gibberella (P < 0.01) (Fig. 4G). The relative abundances of these top three fungi, including Clonostachys, Humicola and Lophiostoma, were also significantly higher in CD (Figs. 4G and 1D). No significant difference was observed in α diversity between these two groups, while PCoA based on weighted UniFrac distance showed the different fungal microflora structures between these two groups (Adonism test, P = 0.04) (Fig. 4H).