In this study, patients in the active stage of UC, according to the Mayo clinical score, were enrolled at the Internal Medicine and Gastroenterology Division at Fondazione Policlinico Universitario “A. Gemelli” IRCCS Hospital.

The inclusion criteria were as follows: diagnosis of UC, mild to moderate active disease (Mayo clinical score, MCS ≤4) or moderate to severe active disease >4, naive to biologic therapies or having failed no more than one line of biologic treatment, and candidate to second-generation therapies. The exclusion criteria were the presence of infectious, ischemic, and actinic colitis or other significant comorbidities.

Healthy subjects, selected for gender and age matching with UC patients, were enrolled at the Human Microbiome Unit of Bambino Gesù Children’s Hospital in Rome, during an epidemiological survey. Subjects with a family history of autoimmune or IBD diseases, with gastrointestinal diseases, and using either antibiotics or pre-/probiotics in the previous 2 months from enrollment were excluded. Fecal samples were collected and stored to generate a reference sample biobank of healthy subjects (BBMRI Human Microbiome Biobank, OPBG).

This study was approved by the Ethics Committee of Fondazione Policlinico Universitario “A. Gemelli” IRCCS Hospital (Protocol 25062019 n.884) and of Bambino Gesù Children’s Hospital, IRCCS (healthy subjects: No. 1113_OPBG_2016), and was conducted in accordance with the Principles of Good Clinical Practice and the Declaration of Helsinki. All participants provided written informed consent for participation in this study.

For bacterial metagenomic analysis, 200 mg of stools were submitted to DNA extraction by QIAmp Fast DNA Stool mini kit (Qiagen, Germany), according to the manufacturer’s instructions.

A 16S rRNA gene fragment, comprising the V3 and V4 hypervariable regions, was amplified using primers reported in the MiSeq rRNA Amplicon Sequencing protocol (Illumina, San Diego, CA, USA).

The approach used for fungal metagenomic analysis started with the lysis step obtained by the incubation of 200 mg of stools resuspended in 500 μl of lysis solution (50 mM of Tris [pH 7.5], 10 mM of EDTA, 28 mM of 2-mercaptoethanol, 10 U/ml of lyticase) (Merck KGaA, Darmstadt, Germany) at 37°C for 30 min in agitation at 850 rpm. After the lysis step, DNA extraction was obtained as described above.

The ITS2 region of approximately 350 bp was amplified using the primers ITS2 5′-GTGARTCATCGAATCTTT-3′ and 5′-GATATGCTTAAGTTCAGCGGGT-3′ (Lemoinne et al., 2020; Del Chierico et al., 2024) under the following conditions: 94°C for 2 min, 35 cycles of 15 s at 94°C, 52°C for 30 s, and 72°C for 45 s. The PCR products were purified using AMPure XP Beads (Beckman Coulters, Brea, CA, USA). A second step of PCR was performed with Illumina-adapted ITS2 primers: 5′-TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGGTGARTCATCGAATCTTT-3′ and 5′-GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGATATGCTTAAGTTCAGCGGGT-3′, following the previously reported PCR conditions applied for 15 cycles.

After a second step of DNA purification, an amplicon PCR indexing step was performed (Nextera XT Index Kit, Illumina) (Del Chierico et al., 2024).

Both bacterial and fungal final libraries were quantified by Quant-iT™ PicoGreen^® dsDNA assay kit (Thermo Fisher Scientific, Waltham, MA, USA), pooled, and sequenced on an Illumina MiSeq™ platform, according to the manufacturer’s specifications. For all amplification steps, negative and positive controls were used to exclude eventual internal and external contaminations for both 16S rRNA and ITS2 sequencing approaches.

All raw sequencing reads are available at the NCBI BioProject database (PRJNA996768 and PRJNA996917) (https://submit.ncbi.nlm.nih.gov/subs/sra/).

Bioinformatics analysis was performed by QIIME2 v.2022.2 software (Bolyen et al., 2019), using DADA2 (Callahan et al., 2016) plugin for quality check, trimming of forward and reverse fastq files, denoising, chimera filtering, and merging reads. The representative sequences of each amplicon sequence variant (ASV) produced with a cutoff of 99% similarity were annotated by using a naive Bayes classifier against the Greengenes reference database (v13.8, http://www.greengenes.secondgenome.com) (DeSantis et al., 2006) for bacteria and the UNITE ITS dynamic database (v9.0, https://unite.ut.ee) (Nilsson et al., 2019) for fungi.

Statistical analyses were performed with R software v4.0.4. For α- and β-diversity analyses, rarefaction was applied on the feature table with absolute frequency, filtering out 16,140 and 2,870 ASVs for 16S and ITS2, respectively. Statistical analyses on α-diversity indices were performed using the non-parametric Mann–Whitney and Kruskal–Wallis tests. The PERMANOVA test was applied to β-diversity matrices. For further analyses, a raw feature table was normalized with the cumulative sum scaling (CSS) method (Paulson et al., 2013).

To characterize and quantify volatile organic compounds (VOCs), 119 stool samples were analyzed with gas chromatography (GC) combined with mass spectrometry (MS) coupled with solid-phase microextraction (SPME) (Douny et al., 2019). The carboxen-polydimethylsiloxane-coated fiber (CAR-PDMS) (85 μm) and the manual solid-phase microextraction (SPME) (Supelco Inc., Bellefonte, PA, USA) were exposed to each sample, for 45 min. The latter was then inserted into the GC injection port for the desorption phase of the samples for 10 min, and GC–MS analyses were carried out using the Agilent Technologies 7890B GC coupled to a 5977A mass selective detector by operating in electron impact mode, equipped with an Agilent DB-HeavyWaX (60 m length, 0.25 mm ID, 0.25 µm) capillary column. All processes were performed under the same conditions reported by Vernocchi et al. (2020).

A match probability of 80%, or greater, was used for VOC identification followed by manual visual inspection of the fragment patterns when necessary. Furthermore, 4-methyl-2-pentanol (final concentration, 400 ppm) was used as an internal standard (IS) in all analyses to quantify the compounds via interpolation of the relative areas in comparison to the IS area (expressed as mg/kg).

VOCs were identified by using a two-step process: the peak spectrum was tested against the NIST (NIST version 2005, NIST 14MS database; National Institute of Standards and Technology, Rockville, MD, USA) mass spectral library database and literature (Garner et al., 2007), and thereafter, in case further confirmation is needed for the NIST identification, it was confirmed by comparing the retention times of the peaks of interest vs. retention times obtained for the reference standards.

ASVs (bacteria: present less than 25% of THE total samples and with relative abundance <0.01; fungi: present less than 25% of the total samples) and metabolites (present less than 10% of the total samples) were filtered out.

Metadata distribution was analyzed by the Shapiro–Wilk test. Gender, age, and clinical features (i.e., corticosteroid therapy, previous therapy, and failure to previous therapies) were evaluated as confounding factors by microbiomeMarker v1.6.0 R package. Age was tested further with the PERMANOVA test by adonis2 function of Vegan v2.6 R package and by linear discriminant analysis (LDA) effect size (LEfSe) (Segata et al., 2011) to exclude it as a confounding factor.

To evaluate the differences in the microbiota, mycobiota, and metabolic profiles between the control (CTRL) and UC patient groups, univariate and multivariate approaches have been applied: linear discriminant analysis effect size, principal component analysis (PCA), and partial least squares-discriminant analysis (PLS-DA). Bacteria function profile was predicted with the Phylogenetic Investigation of Communities by Reconstruction of Unobserved States of Correlation 2 (PICRUSt2) software (Douglas et al., 2020) on 16S rRNA metagenomic data. The p-values were corrected by the Benjamini and Hochberg method to control the false discovery rate (FDR). In all statistical analyses, differences with an adjusted p-value <0.05 were considered significant.

Correlation networks between bacterial and fungal communities were built using Spearman’s correlation by means of the graph and corrr R packages (v1.74.0 and v0.4.4, respectively).

PCA was applied to bacterial, fungal, and metabolite matrices. Then, differential -omic features were screened using variable importance in the projection (VIP) values >1 of the first two principal components of the PLS-DA model and compared with those obtained by differential log fold change of univariate analysis by the Mann–Whitney U test.

The integration of multiple omic data was performed with multivariate approaches: unsupervised ComDim (Common Dimension) multiblock method (Qannari et al., 2000), unsupervised multiblock principal component analysis (MBPCA) (Tchandao Mangamana et al., 2019), and supervised multiblock partial least squares-discriminant analysis (MBPLS-DA) (Brandolini-Bunlon et al., 2019). Each omic matrix (data blocks) was normalized with Frobenius’s method, to harmonize concentration values of metabolites with the relative abundance of microorganisms (Curtasu et al., 2020). Finally, all matrices were joined into a final unique matrix to perform multivariate analyses. MBPLS-DA was validated with the area under the receiver operating characteristic (AUROC) curve, RMSE, Q ², and R ² values.

In this observational study, 53 patients with UC and 37 healthy subjects (CTRL) were enrolled. The cohort characteristics are reported in Table 1 .

By metataxonomic analysis of fecal samples, we obtained 10,729 bacterial and 1,484 fungal ASVs.

We tested gender, age, and some clinical variables as confounding factors in our analyses. No differences in gut microbiota taxa and metabolite distribution were observed for gender, corticosteroid therapy, previous therapy, and failure of previous therapy ( Supplementary Table 1 ). On the contrary, age resulted as a confounding factor for gut microbiota analysis (p-value = 0.032). Then, we performed the beta-diversity analysis on the microbial, fungal, and metabolic variables of subjects grouped by median age (46 years) by the Bray–Curtis algorithm. The PERMANOVA test applied to the dissimilarity matrix showed no statistical significance between the two groups (p-value > 0.05) ( Supplementary Figure 1 ). This result was also confirmed by the LEfSe univariate analysis that revealed no differences in bacterial, fungal, and metabolomic distribution by median age (FDR > 0.05) ( Supplementary Table 2 ).

Comparing UC and CTRL, the α-diversity analysis, calculated on the bacterial composition, revealed no statistically significant differences between UC and CTRL ( Supplementary Figure 2A ), while a significant decrease of the Shannon–Wiener index, calculated on the fungal dataset, was found (p-value = 0.01) ( Supplementary Figure 2B ).

The β-diversity analysis based on the Bray–Curtis algorithm revealed two distinct clusters for UC and CTRL (PERMANOVA p-value < 0.05) ( Supplementary Figures 2C, D ) and an increase of intragroup distance in UC than in CTRL (p-value = 0.0001), for both bacterial and fungal ecosystems ( Supplementary Figures 2E, F ).

To investigate the differences in gut microbiota and mycobiota composition between UC and CTRL, we applied multivariate and univariate approaches ( Figure 1 ).

By the fusion of these results, we assigned Streptococcus, Bifidobacterium, Enterobacteriaceae, TM7-3, Granulicatella, Peptostreptococcus, Lactobacillus, Veillonella, Enterococcus, Peptoniphilus, and Gemellaceae to UC and Akkermansia, Ruminococcaceae, Ruminococcus, Gemmiger, Methanobrevibacter, Oscillospira, Coprococus, and Christensenellaceae to CTRL gut microbiota ( Figures 1A–D ).

For fungi, PCA identified Didymellaceae, Saccharomycetales, Malassezia, Wickerhamomyces, Cutaneotrichosporon, Saccharomyces, Clavispora, Alternaria, and Candida as major fungal markers ( Figure 1E ). PLS-DA analysis revealed the distinctive fecal fungal markers associated with UC ( Figures 1F, G ). Univariate analysis showed a predominance of fungal markers such as Clavispora, Vishniacozyma, and Quambalaria in the CTRL group compared to the UC group, confirming a remarkable difference in the fecal mycobiota between the two groups ( Figure 1H ).

The area under the ROC curve (AUROC) was 0.9393 for bacteria and 0.9951 for fungi, indicating that the applied models have high accuracy in group classifications ( Supplementary Figures 3 , 4 , respectively).

The global composition of the gut microbiota of the study cohorts, represented by the distribution of the ASVs at the phylum, family, and genus levels, is shown in Supplementary Figure 5 . Compared with healthy subjects, the following bacterial markers were more representative in the microbiota of UC patients: Actinobacteria, Proteobacteria, and Bacteroidetes at the phylum level ( Supplementary Figure 5A ); Ruminococcaceae, Enterobacteriaceae, Bifidobacteriaceae, and Streptococcaceae at the family level ( Supplementary Figure 5B ); and Streptococcus, Faecalibacterium, Bifidobacterium, and Bacteroides at the genus level ( Supplementary Figure 5C ). Furthermore, the fungi Ascomycota, Chytridiomycota ( Supplementary Figure 5D ), Saccharomycetaceae, Pleosporaceae, Didymellaceae ( Supplementary Figure 5E ), and Saccharomyces, Malassezia, and Alternaria ( Supplementary Figure 5F ) were the main components of the gut mycobiota in UC patients.

To assess the microbial metabolic pathways, inferred by 16S rRNA sequences, we performed the prediction of pathways of the two cohorts, shown in the LDA plot ( Supplementary Figure 6 ). Twenty pathways, belonging to nine defined metabolic classes and one undefined one, have been associated with the UC profile. Of these, the following pathways were increased in UC: amino acid biosynthesis, aspartate superpathway, carbohydrate degradation, enzyme cofactor biosynthesis, fermentation to lactate, fermentation to lactate/acetate, generation of precursor metabolites and energy, generation of precursor metabolites and energy, purine nucleotide biosynthesis, purine nucleotide de-novo biosynthesis, and terpenoid biosynthesis.

To gain insight into the relationship between taxa from different kingdoms and to gain a more comprehensive understanding of the microbial ecosystems, we applied the interkingdom correlation network analyses to the bacterial and fungal profiles in both cohorts ( Figures 2A, B ).

The UC cohort network was characterized by eight nodes connected by five edges, a relative connectedness of 0.625, and an average number of neighbors of 2.25 ( Figure 2A ). The CTRL cohort network was characterized by 63 nodes connected by 82 edges, a relative connectedness of 1.30, and an average number of neighbors of 3.6 ( Figure 2B ). All significant correlations between bacteria and fungi are listed in Supplementary Table 3 .

We identified 95 filtered molecules by the metabolomic analysis of fecal samples. PCA analysis revealed the presence of two distinct metabolic profiles in UC patients and CTRLs, consisting mainly of 24 metabolic markers ( Figure 3A ).

The PLS-DA analysis showed a low RMSE value (0.329), indicating a high accuracy of the model in predicting the subject classification and highlighting more molecules associated with CTRLs with UC ( Figures 3B, C ). The VIP features are shown in Supplementary Figure 7A , and the AUROC was 0.9944 ( Supplementary Figure 7B ). Finally, the univariate analysis confirmed the results reported by the multivariate approaches ( Figure 3D ). The combination of multivariate and univariate test results highlighted the increase of phenylethyl alcohol and the decrease of cyclopentadecane; 5-octadecene; 5-hepten-2-ol, 6 methyl; 6-methyl-3,5-heptadiene-2-one; hexadecane; 2-butenal; caryophyllene; indole, 3-methyl-; 3-carene; p-Cresol; and 5-hepten-2-one, 6-methyl- in UC.

To reduce the complexity of these multi-omic results, we finally applied three integrative multi-omic approaches to the three omic datasets. The first two were predictive analyses, based on an unsupervised MBPCA and a supervised MBPLS-DA, able to predict discriminant variables (loadings), maintaining multi-omic data separated. The third one was an exploratory approach, carried out by multivariate unsupervised ComDim, in which the three omic matrices were integrated before the analysis. As reported in Figure 4A , the MBPCA identified two distinct gut bacterial, fungal, and metabolic profiles in UC and CTRL.

The MBPLS-DA analyses identified loadings for each sample set ( Figures 4B, C ). Regarding UC, the MBPLS-DA revealed an increase of Bifidobacterium; Streptococcus; TM7-3; 1-hexanol, 2-ethyl-; phenol; benzaldehyde; methyl isobutyl ketone; and 2-heptanone, 4-methyl.

In CTRL, we found the increase of Akkermansia; Gemmiger; Coprococcus; Ruminococcaceae; Christensenellaceae; Clavispora; Vishniacozyma; cyclopentadecane; 3-carene; 1-tridecene; hexadecane; indole, 3-methyl; hexanoic acid ethyl ester; caryophyllene; alpha-pinene; anethole; 2-butenal; p-Cresol; pentadecanal; 2-tetradecene; and 1-heptanol. The RMSE and R ² values were 0.128 and 0.935, respectively, indicating a high performance of this model. The bar plot of VIPs from MBPLS-DA is shown in Figure 4D . The ROC analysis ( Figure 4E ) revealed an AUROC value of 0.9947, indicating a high accuracy of the prediction model.

The application of the ComDim analysis on the three omic matrices integrated confirmed the presence of two distinct UC and CTRL profiles, characterized by bacterial, fungal, and metabolic markers ( Figures 5A, B ).