Graphical abstract graphic file with name fx1.jpg [35]Open in a new tab Highlights * • Gut microbiome shifts in stool precede those in biopsies, indicating faster response * • ESR and SES-CD track bacterial changes, while PLT and BMI reflect fungal variation * • Anaerostipes depletion is linked to the relapse of pCD * • Individual effects exist, with stool samples best representing clinical variations __________________________________________________________________ Pediatrics; Microbiome Introduction The gut microbiota, a diverse and dynamic assemblage of microorganisms residing in the gastrointestinal (GI) tract, plays a critical role in maintaining human health through its intricate interactions with the intestinal epithelium, immune system, and host metabolism.[36]^1 However, disruptions in microbial equilibrium, commonly referred to as dysbiosis, have been implicated in the pathogenesis of various chronic disorders, particularly those affecting the digestive, nervous, and immune systems. Of them, Crohn’s disease (CD) represents one of the most complex inflammatory bowel disease (IBD), characterized by chronic, relapsing inflammation with transmural involvement, capable of affecting any segment of the GI tract.[37]^2 The pathophysiology of CD involves an intricate interplay of genetic susceptibility, environmental triggers, and an aberrant immune response to gut microbiota, making it a lifelong disease requiring sustained therapeutic management.[38]^3 Unresolved intestinal inflammation often leads to irreversible complications, including fibrosis, strictures, fistulas, and malabsorption syndromes, which severely impair quality of life and elevate the risk of adverse outcomes, including malignancies.[39]^4 Recent epidemiological studies indicate an alarming rise in CD incidence among pediatric populations,[40]^5^,[41]^6 reflecting trends observed in adults but with a notably steeper trajectory. This has been further substantiated by Kuenzig et al., who reported increasing global incidence and geographic variability in pediatric-onset CD, emphasizing the need for tailored research in this vulnerable population.[42]^7 Pediatric-onset CD (pCD) is distinct from adult-onset CD in several key aspects, including earlier and more aggressive disease presentation, a higher risk of complications, and the potential for lifelong disease progression. Unlike in adults, where therapy often focuses on symptomatic relief, pediatric patients frequently endure more severe and widespread disease involvement, necessitating a greater reliance on immunomodulatory agents and biologics to achieve fast remission. Moreover, the inflammatory milieu in pCD can disrupt metabolic homeostasis, leading to malnutrition, growth failure, and delayed puberty, adding another layer of clinical complexity. Long-term pharmacotherapy in pCD must be carefully balanced to minimize side effects while effectively controlling inflammation. To improve therapeutic decision-making, it is essential to identify biomarkers that can accurately reflect disease severity, predict relapses, and inform treatment responses. While established clinical biomarkers such as C-reactive protein (CRP) and fecal calprotectin are widely used for disease monitoring,[43]^8 their sensitivity and specificity in predicting long-term outcomes remain suboptimal. Given the heterogeneous nature of pCD, there is an increasing demand for integrative, multi-modal biomarkers that can be reliably applied across diverse patient populations. The gut microbiome has emerged as a key player in the pathophysiology of IBD, with multiple studies demonstrating alterations in bacterial and fungal communities in CD patients,[44]^9^,[45]^10^,[46]^11 including functional shifts linked to inflammation and mucosal integrity following treatments such as exclusive enteral nutrition (EEN) and anti-tumor necrosis factor (TNF) therapy,[47]^12 further supporting its potential role as a biomarker of gut health. However, the choice of sample types introduces significant variability into microbiome studies, limiting the reproducibility and generalizability of findings. A meta-analysis of 13 IBD studies highlighted that variations in sample types is a major confounding factor in disease variability,[48]^13 reinforcing the need for systematic comparisons between different sample types to derive robust conclusions. Despite this, many studies overlook these discrepancies, limiting the ability to establish consistent and reproducible biomarkers across different GI sites. Recent work by Shen et al. further highlighted the need for GI site-specific profiling to better understand microbial heterogeneity in pCD and its clinical relevance.[49]^14 Given these challenges, a comprehensive understanding of the gut microbiome and mycobiome in pCD including variations across different GI sites before and after treatment, is crucial for identifying reliable microbial and clinical markers that can guide precision medicine approaches in pediatric IBD management. In this study, we aim to elucidate the gut microbiota and mycobiota dynamics in pCD by systematically analyzing bacterial and fungal composition across various GI sites, including the colon and terminal ileum, and stool samples from 30 pre-treatment and 27 post-treatment pCD patients, alongside stool samples from 36 healthy individuals. Additionally, we seek to identify clinically relevant biomarkers that correlate with microbial alterations and can be consistently monitored to predict disease progression, remission, and relapse. By providing a comprehensive, site-specific characterization of microbial communities in pCD, this study aspires to enhance disease monitoring strategies, improve biomarker reliability, and ultimately inform personalized treatment approaches for this vulnerable population. Results Microbiome shifts with treatment in stool and biopsy samples A cohort of 30 pediatric patients diagnosed with CD and 36 age-matched healthy controls participated in this study (refer to [50]Table 1). All patients exhibited moderate to severe CD, as quantified by the PCDAI. To identify disease-specific traits, diagnostic, and prognostic markers that are consistent across various sample groups, we collected biological specimens from the stool, colon, and terminal ileum. The study comprised 30 pediatric patients and age-matched 36 healthy controls (see Methods). Additionally, samples were collected from 27 subjects after treatment. Samples from the healthy control stool, patient stool, patient colon tissue, and patient terminal ileum tissue were designated as “healthy control”, “patient St”, “patient C″, and “patient Ti”, respectively. Comparisons of clinical indices between healthy controls and pre-treatment samples revealed significantly lower levels of inflammatory markers in the healthy group across all indices. Similarly, comparisons between pre-treatment and post-treatment samples showed a significant reduction in inflammation following treatment (p < 0.05 for clinical, endoscopic index and laboratory markers). While post-treatment samples exhibited no significant differences in CRP and WBC levels compared to healthy controls (CRP: p = 0.191, WBC: p = 0.687), ESR and fecal calprotectin levels remained significantly elevated in the post-treatment group (ESR: p = 0.007, Fecal calprotectin: p = 0.014). To examine the bacterial microbiome and fungal mycobiome across these sample types and subject groups, amplicon sequencing was successfully performed and achieved an average of 62,116 (16S) and 72,632 (ITS2) high-quality microbial reads ([51]Figure S1). The resulting number of reads after each step was also shown ([52]Table S2). Furthermore, an array of clinical data including BMI, PCDAI, SES-CD, ESR, CRP levels, fecal calprotectin levels, WBC, and PLT was documented. Table 1. Demographic and clinical characteristics of 30 pediatric Crohn’s disease patients and healthy controls Subjects Healthy control (n = 36) Pre-treatment (n = 30) Post-treatment (n = 27) p value (HC vs. Pre-) p value (HC vs. Post-) p value (Pre- vs. Post-) Age (yr) [Median IQR] 14.51 [13.2–15.8] 13 [12–15.5] 14 [13–15.5] 0.887 0.726 0.617 __________________________________________________________________ Disease location __________________________________________________________________ L1 [Terminal ileal/limited cecal] (n, %) – 4 (13.3%) 1 (3.7%) – – 0.206 L2 [Colonic] (n, %) – 1 (3%) 3 (11.1%) – – 0.259 L3 [Ileocolonic] (n, %) – 25 (83.3%) 3 (11.1%) – – <0.001 Upper GI involve (n, %) 24 (80%) 6 (22.2%) – – <0.001 __________________________________________________________________ Gender __________________________________________________________________ Male (n, %) 27 (75%) 24 (75%) 21 (77.8%) 0.629 0.798 0.837 Female (n, %) 9 (25%) 6 (25%) 6 (22.2%) – – – __________________________________________________________________ Anthropometric data __________________________________________________________________ Weight (kg) [Median IQR] – 48.2 [37–61] 56 [48.2–64.9] – Height (cm) [Median IQR] – 159 [150–170] 165.8 [160.25–171.6] – BMI (kg/m^2) [Median IQR] 21.1 [18.9–24.1] 18.4 [16.1–20.4] 20 [18.5–22.15] 0.016 0.199 0.153 __________________________________________________________________ Clinical index __________________________________________________________________ PCDAI [Median IQR] – 30 [25–35] 2.5 [0–5] – – <0.001 __________________________________________________________________ Endoscopic index __________________________________________________________________ SES-CD [Median IQR] – 20 [14.3–27] 0 [0–1.5] – – <0.001 __________________________________________________________________ Laboratory markers __________________________________________________________________ WBC [Median IQR] 6090 [5530–7400] 9000 [7502–10620] 6180 [5235–7045] 0.038 0.687 0.024 PLT [Median IQR] – 366 [297–458] 256 [219–283] – – <0.001 ESR [Median IQR] 3 [2–4.75] 28 [18–72.3] 6 [2–13] <0.001 0.007 <0.001 CRP [Median IQR] 0.03 [0.03–0.06] 0.87 [0.32–1.99] 0.06 [0.06–0.115] <0.001 0.191 <0.001 Fecal calprotectin (>1000 μg/g) (n, %) [Median IQR] 38.5 [14.23–62.25] 23 (77%) 144 [32.3–597.5] <0.001 0.014 <0.001 __________________________________________________________________ Treatment __________________________________________________________________ Mesalazine – – 27 (100%) – – – Immunomodulator – – 27 (100%) – – – Exclusive enteral nutrition (EEN) – – 20 (71.4%) – – – Infliximab (IFX) – – 23 (85.2%) – – – Antibiotics within 2 months (n, %) – 1 (3%) 1 (3%) – – – Other drug treatment (n, %) 0 (0%) 0 (0%) 0 (0%) – – – [53]Open in a new tab The cohort includes 30 pCD patients (aged 10–18 years), assessed at pre-treatment (n = 30) and post-treatment (n = 27), alongside 36 Healthy Controls. Data are presented as medians with interquartile ranges (Q1–Q3). Significant differences were observed between pre- and post-treatment groups in inflammatory markers (WBC, PLT, ESR, CRP, fecal calprotectin) and disease activity indices (PCDAI, SES-CD), indicating clinical remission. All patients received mesalazine and immunomodulators; the majority also received infliximab (85.2%) and exclusive enteral nutrition (71.4%). No antibiotics were administered during the sample collection period. p-values reflect group comparisons; unpaired t-tests were used for continuous variables, and chi-square tests were applied to categorical variables (e.g., gender, pre-treatment fecal calprotectin >1000 μg/g), performed using GraphPad Prism v10.0.2. We began by comparing Shannon and beta-diversity across sample groups to obtain an overview of microbial composition ([54]Figures 1A‒1C). In stool samples, the healthy control group had a higher Shannon index compared to both the pre-treatment and post-treatment patient St groups (adjusted p = 0.047, unpaired t test with FDR correction). Interestingly, the pre-treatment patient St group showed greater variance in alpha diversity compared to the post-treatment patient St group (p < 0.001, F-test). Additionally, the average within-group weighted UniFrac distance of patient St samples was higher than that of other groups (healthy control: 0.27, patient St: 0.44, patient C: 0.29, patient Ti: 0.28), indicating greater intra-group variability in patient stool samples ([55]Figure 1B). A significant difference was observed between patient and healthy control groups within stool samples (R^2 = 0.056, p = 0.005, PERMANOVA test). In contrast, biopsy samples showed a different pattern from stool samples ([56]Figures 1B and [57]S2A). Specifically, the Shannon indices for biopsy samples significantly decreased after treatment (adjusted p = 0.0002, unpaired t test with FDR correction), whereas stool samples showed no statistically significant change. Although differences in alpha diversity between pre- and post-treatment samples were modest in visual appearance, beta diversity analysis revealed statistically significant compositional differences across all sample types, including stool, colon, and Ti (p < 0.05, PERMANOVA tests; [58]Table S5). However, the corresponding R^2 values were relatively low (0.043–0.060), indicating that the magnitude of microbiome shift after treatment was subtle but consistent ([59]Figures 1A and [60]S2). Despite locational differences, colon and Ti samples demonstrated similar alpha- and beta-diversity (p = 0.754, unpaired t test; R^2 = 0.004, p = 0.98, PERMANOVA test). Notably, the microbiome of post-treatment patient St samples resembled that of healthy control samples (R^2 = 0.022, p = 0.196), suggesting a recovery toward a healthier microbiome state ([61]Figure 1C). Figure 1. [62]Figure 1 [63]Open in a new tab Characteristics of microbiota in patients and Healthy Controls (A) Boxplots showing alpha-diversity at the ASV feature level, measured by Shannon index for each sample group. Differences between groups were assessed using the unpaired t test with FDR correction (ns: p > 0.05, ∗: p < 0.05, ∗∗: p < 0.01, ∗∗∗: p < 0.001). Beta-diversity plots based on the weighted UniFrac distance metric (B) for healthy controls and pre-treatment groups (Healthy Control vs Patient St; R^2 = 0.0056, p = 0.005, PERMANOVA test), (C) comparison among healthy controls, pre-treatment, and post-treatment groups in fecal microbiome (pre-treatment vs post-treatment; R^2 = 0.043, p = 0.034, PERMANOVA test). Venn diagram depicting the number of core microbiomes at the genus-level, for (D) pre-treatment and (E) post-treatment groups. Average relative abundance of microbiome at the (F) phylum and (G) genus levels, including core phyla and genera present in ≥75% each sample groups, and ≥1% mean relative abundance. Rest of genera were aggregated into the “Others”. To further investigate compositional differences between sample groups, we identified dominant phyla and genera, defining core phylum and genus as those with at least 1% mean relative abundance present in a minimum of 75% of the samples ([64]Figures 1D and 1E; [65]Table S3). Two genera, Bacteroides (prevalence = 92.8%) and Faecalibacterium (prevalence = 68%) were identified as core genera across all pre-treatment sample groups. Anaerostipes and Lachnospiraceae UCG-004 were unique core members in healthy control samples, while Escherichia-Shigella was predominant in patient St samples. The patient St group had the fewest core genera, likely due to reduced alpha-diversity and increased sample variability, as also reflected in the principal coordinates analysis (PCoA) plot ([66]Figure 1B). In contrast, pre-treatment colon and Ti samples shared up to 27 core genera, indicating high within-group similarity ([67]Figure 1D), while in post-treatment samples, the number of core genera greatly decreased ([68]Figure 1E). To compare detailed compositional differences between groups, we calculated the mean relative abundances of each core phylum and genus ([69]Figures 1F and 1G). Proteobacteria, often associated with pathogenic bacteria, was dominant in both stool and biopsy samples from patients (6.0%),[70]^15 while Actinobacteria, crucial for gut homeostasis,[71]^16 was predominant in healthy controls (7.6%). Patient samples including both stool and biopsy, exhibited higher mean relative abundance of Escherichia-Shigella (patient St: 12.2%, patient C: 16.4%, patient Ti: 15.0%) and lower levels of Bifidobacterium (healthy control: 7.4%, patient C: 1.7%, patient Ti: 1.7%) and Faecalibacterium (healthy control: 10.0%, patient St: 4.2%, patient C: 5.4%, patient Ti: 5.4%) in pre-treatment samples. Post-treatment stool samples showed a decrease in Proteobacteria and an increase in Firmicutes. While stool samples showed significant compositional differences at both phylum and genus levels between pre- and post-treatment phases, biopsy samples remained similar at the phylum level ([72]Figure S3). Also, in post-treatment biopsy samples, a considerable abundance of Escherichia-Shigella and Ruminococcus gnavus group was observed, which are frequently associated with CD pathology.[73]^17^,[74]^18^,[75]^19^,[76]^20 Overall, stool and biopsy sample showed different patterns of variation, which underscores the complexity of microbial composition changes in response to treatment. Microbial signatures identify potential biomarkers across sample types To further explore and statistically validate group-specific significant microbial features, we performed a series of differential abundance analyses using linear discriminant analysis effect size (LEfSe), applying a threshold of p value <0.05 and LDA score >3.5. The analysis demonstrated that healthy control samples had a higher abundance of the phyla Firmicutes and Actinobacteria ([77]Figure 2A). At the genus level, these samples were enriched in Agathobacter, Bifidobacterium, Blautia, Dialister, Faecalibacterium, Lachnospira, Roseburia, and Subdoligranulum. We additionally applied microbiome multivariable associations with linear models (MaAsLin2) to validate the results of LEfSe. Here, only Klebsiella and Lachnoclostridium were excluded from enriched genera in healthy control, since MaAslin2 did not identify significant associations of those genera. Conversely, patient samples were characterized by an increased abundance of the phylum Proteobacteria and genera such as Enterococcus, Escherichia-Shigella, and the Ruminococcus gnavus group, which are often associated with inflammation and dysbiosis.[78]^20^,[79]^21^,[80]^22 Figure 2. [81]Figure 2 [82]Open in a new tab Differential abundance analysis using linear discriminant analysis effect size (LEfSe) and relative abundance boxplots of microbial features (A and B) Circular cladogram illustrating microbial species with significant differences among groups. Nodes colored green, red, blue represent key microbial features in healthy control, pre- and post-treatment patient stool samples, respectively. Microbial features with absolute LDA scores greater than 3.5 are shown. (C) Boxplots of relative abundance for each genus showing significant differences between healthy controls and pre-treatment patient stool samples. Upon comparing post-treatment and pre-treatment stool samples, we observed an enrichment of genera such as Anaerostipes, Blautia, Butyricicoccus, Lachnoclostidium, and Tuzzerella, while Enterococcus and Prevotella_9 were depleted ([83]Figure 2B). However, when comparing post-treatment patient stool samples to healthy control samples, fewer distinct microbial features were identified ([84]Figure S4A), indicating a microbiome composition shift toward a healthier state. Our differential abundance analysis revealed no significant differences in microbial composition between colon and Ti samples in either the pre-treatment or post-treatment groups, indicating their similar nature. Interestingly, most genera depleted in patient St samples were not similarly depleted in biopsy samples ([85]Figure 2C), prompting us to compare patient stool and biopsy samples to identify type-specific characteristics. Both pre-treatment and post-treatment biopsy samples were enriched in Escherichia-Shigella and Prevotella compared to stool samples ([86]Figures S2B and S2C). Additionally, pre-treatment biopsy samples compared to stool samples showed enrichment in genera such as Blautia, Faecalibacterium, Klebsiella, and Lachnospira, which were depleted in patient St samples ([87]Figure S4C). Notably, genera typically reported as depleted in CD, such as Akkermansia, Blautia, Faecalibacterium, Lachnospira, and Prevotella were enriched in pre-treatment biopsy samples.[88]^21^,[89]^22 Conversely, genera such as Fusobacterium, Escherichia-Shigella, and Ruminococcus gnavus group often considered pathogenic,[90]^20^,[91]^22^,[92]^23^,[93]^24 were enriched in post-treatment biopsy samples ([94]Figures S4B‒S4E). This observation, together with the decreased alpha diversity of post treatment biopsy samples, suggests delayed microbial shifts in biopsy samples compared to stool samples. To further elucidate functional differences between sample groups, we examined the enriched and depleted pathways in patient St samples using PICRUSt2 and ALDEx2 ([95]Figure S5). Patient St samples were enriched in pathways related to ubiquinol biosynthesis, typically associated with oxidative stress response.[96]^25 In contrast, healthy control St samples were enriched in pathways for multiple biosynthesis (polyamine, L-arginine, phosphatidylglycerol, L-lysine, L-methionine) and degradation, which are linked to anti-inflammatory effects ([97]Figure S5A).[98]^26^,[99]^27^,[100]^28 Additionally, when comparing pre-treatment and post-treatment Patient St samples, pathways such as “galactose degradation” and “methanogenesis from acetate” were enriched in the post-treatment group ([101]Figure S5B). Additionally, to assess the impact of treatment type (e.g., infliximab therapy (IFX) or EEN) on clinical outcomes and microbiome composition after treatment, unpaired t-tests were performed on the delta values of clinical indices between pre- and post-treatment samples in 27 subjects, stratified by treatment regimen. No significant differences were observed in any clinical parameter or Shannon diversity across treatment types (p > 0.05 for all comparisons, [102]Table S4). Similarly, PERMANOVA analysis using beta diversity showed no significant separation by treatment (EEN: R^2 = 0.011, p = 0.498; IFX: R^2 = 0.027, p = 0.056). Alpha diversity also did not differ significantly between treatment groups across any sample type (unpaired t test with FDR correction; adjusted p > 0.05 for all tests, [103]Figure S6). In summary, treatment regimen was not a major determinant of either microbiome composition or clinical index changes during the treatment period. Aspergillus prevalence and fungal variation in samples versus biopsy samples Next, we analyzed the fungal mycobiome across different sample types to elucidate distinct patterns and differences. This analysis involved amplifying ITS2 regions and utilizing the DADA2 pipeline, similar to that used for bacterial microbiome analysis (see Methods). Unlike 16S samples, ITS2 samples exhibited a lower species richness, with a maximum of 20 species detected ([104]Figure S1B). Also, stool samples demonstrated a higher Shannon diversity index compared to colon samples (adjusted p = 0.026; [105]Figure 3A). Significant differences were observed only between stool and colon samples through PERMANOVA results (patient St vs. patient C: R^2 = 0.086, p = 0.007; patient St vs. patient Ti: R^2 = 0.029, p = 0.168; patient C vs. patient Ti: R^2 = 0.029, p = 0.190; [106]Figure 3B). In contrast to bacterial samples, stool and terminal samples didn’t show significant difference according to both alpha and beta diversity analysis. Regarding prevalent genera, Aspergillus (prevalence = 42.17%), Penicillium (prevalence = 30.12%), Cladosporium (prevalence = 28.92%), and Naganishia (prevalence = 27.71%) were prevalent across all sample types. Contrary to 16S microbiome analysis results, patient St samples exhibited a notable overlap in prevalent genera ([107]Figure 3C; [108]Table S3). The fungal mycobiome predominantly comprised two phyla, Ascomycota and Basidiomycota ([109]Figures 3D and [110]S7). Notably, patient Ti samples were uniquely characterized by the presence of Malessezia (Patient Ti; prevalence = 29.6%, mean relative abundance = 8.2%), while patient C samples were distinguished by a presence of Naganishia (patient C; prevalence = 33.3%, mean relative abundance = 11.7%), suggesting specific fungal niches in the gut ([111]Figure 3E). Further analysis using LEfSe and MaAsLin2 identified that only Saccharomyces was significantly enriched in stool samples compared to biopsy samples ([112]Figure S8). Figure 3. [113]Figure 3 [114]Open in a new tab Diversity and composition of the fungal mycobiome in patient samples (A) Boxplots showing the Shannon index for each sample type at the ASV feature level. The unpaired t test with FDR correction was performed (∗: p < 0.05). (B) Beta-diversity plot using the weighted UniFrac distance matrix (Patient St vs Patient C; R^2 = 0.086, p = 0.007, PERMANOVA test). (C) Venn diagram of mycobiome at the genus level with a prevalence threshold of 0.1. Average relative abundance of the mycobiome at the (D) phylum and (E) genus levels. Individual microbial profiles show intra-patient similarities In both the bacterial microbiome and fungal mycobiome, we observed distinct differences between stool and biopsy samples. Procrustes analysis was performed to evaluate the correlation between weighted UniFrac distance matrices of different groups and to statistically confirm individual effects ([115]Figure 4A). Overall, significant correlations were identified between patient stool and biopsy samples in both pre-treatment (St-C: r = 0.58, St-Ti: r = 0.65, C-Ti: r = 0.92) and post-treatment samples (St-C: r = 0.59, St-Ti: r = 0.53, C-Ti: r = 0.87), indicating that different sample types from the same individual exhibit similar microbial profiles ([116]Figure 4B). This finding suggests the presence of individual-specific microbial characteristics across sample types. Interestingly, pre-treatment stool samples demonstrated relatively consistent relationships with all other sample types, including post-treatment samples and clinical data, underscoring their representational strength. Also, compared to biopsy samples, post-treatment stool samples showed slightly higher correlation values (St: 0.48, C: 0.45, Ti: 0.44). In contrast, for fungal samples, only patient colon biopsy (patient C) samples exhibited notable similarities with pre-treatment stool bacterial samples and post-treatment colon bacterial samples. The correlations between bacterial (16S rRNA) and fungal (ITS2) microbiomes were generally weak, suggesting a limited interdependence between bacterial and fungal community dynamics. Figure 4. [117]Figure 4 [118]Open in a new tab Similarity comparisons between sample types through Procrustes analysis based on beta diversity distance (A) Description of Procrustes analysis. (B) Procrustes analysis result based on the weighted UniFrac distance of microbiome and mycobiome. The bubble size reflects the Procrustes correlation coefficient value and color transparency indicates the p value. ESR and SES-CD track bacterial shifts, while PLT and BMI reflect fungal changes. We then examined the relationship between clinical factors and microbial communities across different sample types to identify key markers of microbial variation. Using PERMANOVA and Pearson correlation analyses, we assessed how various clinical factors correlate with microbial states ([119]Figure 5). ESR demonstrated strong correlations with microbial variations in stool, colon, and Ti samples, showing consistently high R^2 values and low p values (St: R^2 = 0.13, p = 0.005; C: R^2 = 0.086, p = 0.043; Ti: R^2 = 0.102, p = 0.019; [120]Table S6). In contrast, the log-transformed PLT was significant only in stool samples (R^2 = 0.135, p = 0.004) but not in biopsy samples ([121]Figure 5A). In post-treatment samples, the SES-CD showed the highest R^2 values and low p values across sample types (St: R^2 = 0.075, p = 0.088; C: R^2 = 0.118, p = 0.028; Ti: R^2 = 0.057, p = 0.152; [122]Figure S9A). However, ITS2 mycobiome analyses yielded inconsistent results ([123]Figure 5B). For beta-diversity, BMI (St: R^2 = 0.097, p = 0.042, C: R^2 = 0.248, p = 0.001, Ti: R^2 = 0.015, p = 0.729), log-scaled PLT (St: R^2 = 0.043, p = 0.264, C: R^2 = 0.111, p = 0.044, Ti: R^2 = 0.021, p = 0.581) and log-scaled WBC (St: R^2 = 0.046, p = 0.272, C: R^2 = 0.009, p = 0.919, Ti: R^2 = 0.097, p = 0.066) accounted for variability in the ITS2 mycobiome. Figure 5. [124]Figure 5 [125]Open in a new tab Associations between clinical factors and gut microbiome composition The significance of each clinical variable on beta-diversity metrics shows (A) for the weighted UniFrac distance for the 16S microbiome and (B) the ITS2 mycobiome. These values are visualized as a bubble heatmap, where bubble size reflects the proportion of variance explained by each clinical variable (R^2), and color transparency denotes the p value. Pearson correlation heatmap showing the relationship between clinical indices and alpha diversity values for 16S microbiome and mycobiome within pre-treatment patients from (C) stool, (D) colon, and (E) terminal ileum samples. To further validate the results of PERMANOVA based on weighted UniFrac distances, we additionally conducted a distance-based redundancy analysis (db-RDA) using Bray-Curtis distance matrices ([126]Table S7). This multivariate approach largely supported our initial findings. In the ITS2 dataset, BMI showed significant associations with fungal community composition in both stool and colon (St: p = 0.046; C: p = 0.026) samples, but not in Ti (p = 0.867). For the post-treatment 16S microbiome, SES-CD was significantly or near-significantly associated with community variation in the colon and terminal ileum (St: p = 0.362; C: p = 0.062; Ti: p = 0.038). In pre-treatment samples, ESR exhibited marginal associations (St: p = 0.429; C: p = 0.088; Ti: p = 0.082), aligning with trends observed in PERMANOVA. Pearson correlation analysis on alpha-diversity indices revealed a negative correlation with ESR in stool, colon, and Ti samples ([127]Figures 5C‒5E). In post-treatment samples, SES-CD, ESR, and CRP consistently exhibited negative correlations with alpha-diversity indices across sample types ([128]Figures S9B‒S9D). Alpha-diversity of the fungal mycobiome showed significant correlations exclusively with BMI but was not consistent across all sample types. Notably, the alpha-diversity of the ITS2 fungal mycobiome negatively correlated with the alpha-diversity of 16S microbiome only in biopsy samples ([129]Figures 5D and 5E), which may imply a competitive relationship between bacteria and fungi in the local gut environment. Additionally, principal component analysis (PCA) of clinical data identified ESR and log-scaled PLT as significant factors in pre-treatment subjects, contributing 19.6% and 19.4% to the first PC1 (32.4%, variance proportion), respectively ([130]Figures S10A and S10B). In post-treatment subjects, ESR, SES-CD, and CRP were the most significant, contributing to the first PC1 (30.7%, variance proportion) ([131]Figures S10C and S10D). However, though BMI was shown to track the variation of the mycobiome, BMI contributed less than 12% to PC1 and PC2, indicating that the BMI was not effectively representing the clinical data. Overall, ESR, SES-CD, and log-scaled PLT as particularly effective markers reflecting variance in clinical data and microbiome/mycobiome. Microbial correlations with clinical indices suggest biomarkers We further examined phylum-level ratios to evaluate their potential associations with clinical indicators, focusing particularly on the Firmicutes/Bacteroidetes (F/B) ratio, which has previously been reported to be reduced in CD ([132]Table S8).[133]^29 In our dataset, the Proteobacteria/Bacteroides ratio and the combined relative abundance of Firmicutes and Bacteroides demonstrated relatively consistent correlations with clinical indices across both pre- and post-treatment samples. In contrast, the F/B ratio itself did not show any significant associations with clinical markers in the pre-treatment samples. Interestingly, however, in post-treatment samples, the F/B ratio exhibited a strong positive correlation with key inflammatory markers, including ESR, SES-CD, and CRP (mean r; St: 0.60, C: 0.38, Ti: −0.01). These findings suggest that biomarker relevance may shift following treatment, and that the F/B ratio may have limited utility for disease monitoring in the post-treatment phase due to its altered correlation pattern. We evaluated the relationships between microbial genera and selected clinical indices, i.e., ESR, SES-CD, log-scaled PLT and BMI to identify biomarkers for disease severity. Using MaAsLin2, we generated mixed effect models to calculate correlations of genera with these clinical indices, visualizing the results with a heatmap in a Circos plot ([134]Figure 6). Genera from pre-treatment samples showed more correlations with clinical indices compared to those from post-treatment samples, likely due to a significant reduction in clinical indices following treatment ([135]Table 1). In pre-treatment samples, genera Akkermansia, Lachnoclostridium, Lachnospira, Phascolarctobacterium, Roseburia, Ruminococcus gnavus group, and Subdoligranulum exhibited negative correlations with ESR while Fusobacterium had positive correlation with ESR in more than two different sample types. In post-treatment samples, Anaerostipes, Blautia, and Faecalibacterium were negatively correlated with ESR or SES-CD across different sample groups. Overall, genera from the phylum Bacteroidota and family Lachnospiraceae were negatively correlated with clinical indices. Also, Enterococcus and Klebsiella were positively correlated with ESR and SES-CD in pre-treatment patient stool samples. In the fungal mycobiome, Malassezia represented positive correlations with ESR and log-scaled PLT and negative correlations with BMI, while Naganishia and Cladosporium showed negative correlation with BMI in stool samples, indicating their relevance to the severity of pediatric CD ([136]Figure 6). Figure 6. [137]Figure 6 [138]Open in a new tab Circos plot illustrates inter- and intra-correlation of genera, including a circular heatmap of correlation coefficients from MaAsLin2 analysis The outermost layer 3 and the subsequent layer 2 represent the phylum- and family-level of the genera, respectively. The circular heatmap (layer 2) displays MaAsLin2-calculated correlation coefficients for ESR, SES-CD, log-scaled PLT, and BMI (from outer to inner layers) in relation to bacterial and fungal genera. The inner layer 1 lists genera grouped by sample type, with color coding shown at the bottom of the panel. The innermost network layer depicts significant correlations between genera, as determined by “SpiecEasi” (p < 0.05). Red and blue lines indicate positive and negative correlations, respectively, with lines connecting the same genus excluded from the plot. To further investigate individual correlations among microbial and fungal genera, we analyzed co-occurrence networks, identifying keystone genera within each sample group. In pre-treatment samples, the network comprised 58 genera, with Bacteroides, Enterococcus, Lachnoclostridium, and Veillonella emerging as key nodes, each demonstrating connections with more than four edges. In contrast, the network from post-treatment samples included 88 genera, with Faecalibacterium, Helicobacter, Neisseria, and Ruminococcus gnavus group identified as prominent nodes. The shift in keystone genera between pre-treatment and post-treatment samples highlights the dynamic nature of the microbiome in response to therapeutic interventions. Anaerostipes depletion is linked to disease relapse Lastly, we compared the post-treatment samples to assess differences based on disease relapse among subjects. Out of 27 patients, 5 experienced a relapse. LEfSe and MaAsLin2 analysis revealed a consistent depletion of the genus Anaerostipes across all sample groups in patients who relapsed ([139]Figures S11A‒S11C). A PERMANOVA test revealed a significant difference in beta-diversity between relapsed and non-relapsed groups (R^2 = 0.036, p = 0.02), with colon samples showing particularly pronounced differences (R^2 = 0.095, p = 0.04) ([140]Figure 7A). Interestingly, comparing the relative abundance of Anaerostipes between pre- and post-treatment relapse samples revealed a notable deepening of its depletion over time (pre-treatment, St: adjusted p = 0.028, C: ns, Ti: ns; post-treatment, St: adjusted p = 0.028, C: adjusted p = 0.028, Ti: adjusted p = 0.041) ([141]Figure 7B). This suggests that Anaerostipes depletion may be associated with disease relapse and could potentially serve as a biomarker for monitoring CD progression. Figure 7. [142]Figure 7 [143]Open in a new tab Anaerostipes depletion is associated with relapse in pediatric Crohn’s disease (A) Beta-diversity based on weighted UniFrac distances of post-treatment microbiome samples, showing distinct clustering between patients with clinical relapse (red) and those in remission (green) (Relapse vs Remission; R^2 = 0.036, p = 0.020, PERMANOVA test). (B) Relative abundance of the Anaerostipes across sample types and treatment stages. A significant depletion of Anaerostipes was observed in patients who experienced relapse, particularly in post-treatment stool and biopsy samples. Asterisks denote statistical significance (∗: p < 0.05, ns: not significant; Welch’s t test with FDR correction). Discussion This study provides a comprehensive multi-compartmental analysis of bacterial and fungal communities in pCD, revealing how sample type and treatment status shape microbial patterns and their associations with clinical indicators. By profiling the microbiome and mycobiome across stool, colon, and terminal ileum (Ti) samples before and after treatment, we identified spatial and temporal heterogeneity in microbial responses, reinforcing the critical role of sample selection in microbiome-based biomarker discovery. Stool and mucosal samples displayed distinct ecological signatures, with biopsy samples showing delayed microbial shifts post-treatment. This lag may reflect the increased resilience of mucosa-associated communities due to their proximity to the host epithelium and immune cells, and their adaptation to localized conditions such as oxygen gradients and nutrient availability.[144]^30^,[145]^31^,[146]^32 These findings align with previous studies showing that stool samples, as composites of luminal and mucosal microbiota, more rapidly reflect treatment-induced changes, and may offer a practical advantage for non-invasive monitoring of disease activity.[147]^33 Taxonomic analysis further revealed that taxonomic resolution and sample type are critical when interpreting treatment response. For example, while alpha diversity remained largely unchanged in stool samples after treatment, phylum-level shifts suggested a complex compensatory dynamic. Firmicutes-associated genera increased while richness within Bacteroidota and Proteobacteria decreased ([148]Table S9), resulting in a stable but compositionally distinct community. Such phylum-level shifts have been reported previously,[149]^34^,[150]^35 suggesting that newly dominant amplicon sequence variants (ASVs) in Firmicutes were insufficient to notably enhance diversity. Despite unchanged relative abundance of Bacteroidota (p = 0.929), ASVs richness within this phylum decreased significantly (p = 0.002), contributing to reduced overall alpha diversity. At the genus level, the depletion of Prevotella and Prevotella_9, known for both pro-,[151]^36 and anti-inflammatory[152]^37 properties depending on species context, emphasizes the limitations of genus-level interpretation and the need for higher taxonomic resolution to distinguish strain-specific effects. Consistent with previous findings,[153]^22 our analyses of microbial diversity, composition, and differential abundance revealed clear distinctions between stool and biopsy samples. In particular, patient St samples showed reduced alpha diversity and marked intra-individual variability,[154]^38^,[155]^39^,[156]^40 while post-treatment microbial communities began to resemble those of healthy controls. However, the absence of significant alpha diversity improvement in stool samples contrasts with earlier reports of post-treatment recovery,[157]^41^,[158]^42 possibly reflecting the aforementioned taxonomic compensations. LEfSe analysis identified enrichment of Enterococcus and Escherichia-Shigella in pCD patients, while short chain fatty acid (SCFA)-producing genera such as Lachnospira, Roseburia, and Subdoligranulum were depleted.[159]^43^,[160]^44^,[161]^45^,[162]^46^,[163]^47 Notably, six genera (i.e., Agathobacter, Bifidobacterium, Dialister, Roseburia, and Subdoligranulum) were consistently depleted in patients and have been implicated in previous pCD studies,[164]^48 further supporting their potential as biomarker. Meanwhile, the observed enrichment of Ruminococcus gnavus group, which is a genus frequently linked to IBD or CD,[165]^49^,[166]^50^,[167]^51 contrasted with its negative correlation with ESR in biopsy samples ([168]Figure 6), emphasizing the context-dependent behavior of even well-established taxa. At the functional level, our analysis of predicted metabolic pathways suggested microbial adaptation to host inflammation, with ubiquinol biosynthesis enriched in pCD patients, possibly reflecting microbial oxidative stress responses.[169]^51 In contrast, pathways associated with anti-inflammatory metabolites (e.g., L-arginine, L-lysine, L-methionine, phosphatidylglycerol, and polyamine) were predominant in healthy controls, corroborating their reported immunomodulatory roles ([170]Figure S5A).[171]^26^,[172]^27^,[173]^28 The enrichment of “methanogenesis from acetate” in post-treatment stool samples may also indicate microbial contributions to improved gut environment through modulation of luminal gas pressures ([174]Figure S5B).[175]^52 In the mycobiome, we found that site-specific patterns, with dominant taxa in biopsy samples (e.g., Aspergillus, Naganishia, Penicillium, and Cladosporium) differing from those in stool ([176]Figure 3C). The relative paucity of Saccharomyces and Candida in biopsy specimens, despite their frequent identification in gut mycobiome studies,[177]^53^,[178]^54 may be attributed to heightened mucosal inflammation. Importantly, stool-derived mycobiomes exhibited stronger correlations with clinical indices, suggesting their utility in non-invasive disease monitoring. Mycobiome diversity was lower than bacteriome diversity, limiting overlap across samples but reinforcing the existence of niche-specific fungal communities. Among clinical variables, ESR and SES-CD were most strongly associated with bacterial microbiome variation across sample types ([179]Figure 6). While fecal calprotectin is a widely considered a sensitive marker for IBD,[180]^55 its correlation was weaker in our pre-treatment dataset, potentially due to its sporadic availability. ESR’s strong correlation and stability may be attributable to its integration of long-term inflammatory burden, less susceptible to short-term perturbations compared to CRP or calprotectin ([181]Figure 5).[182]^56^,[183]^57^,[184]^58^,[185]^59 Moreover, SES-CD, an endoscopic measure of mucosal healing,[186]^60 was especially effective in reflecting microbial shift in biopsy samples. For the mycobiome, log-scaled PLT, and BMI emerged as moderate correlates, consistent with prior studies linking fungal composition to host nutritional status and inflammatory load.[187]^61^,[188]^62^,[189]^63^,[190]^64^,[191]^65 Fungal taxa have been implicated in modulating or responding to host metabolism through mechanisms such as fermentation of dietary polysaccharides, production of SCFAs, and regulation of intestinal barrier function, processes tightly linked to nutritional homeostasis. Therefore, among various clinical indices, BMI, being a direct indicator of nutritional and metabolic state,[192]^66 may plausibly exhibit the strongest association with fungal composition. Interestingly, Lachnoclostridium and Ruminococcus gnavus group demonstrated divergent associations with clinical indices.[193]^22^,[194]^23 The depletion of Lachnoclostridium, contrary to prior IBD reports[195]^23^,[196]^24 may reflect its context-specific immunomodulatory role as suggested by recent gene cluster analyses.[197]^67 Ruminococcus gnavus group demonstrated a negative correlation with ESR in biopsy samples but was enriched in pre-treatment stool samples and all post-treatment samples. This discrepancy could be attributed to the delayed reflection of disease conditions in biopsy samples, particularly given its enrichment in post-treatment samples. However, the substantial differences in genera associated with clinical indices before and after treatment suggest that different biomarkers may be required depending on the disease stage and patient condition. Despite the challenges in identifying fungal biomarkers due to the lack of healthy control samples, we found that Cladosporium and Malassezia, which are known to be associated with CD,[198]^68 underscoring their significance in pCD. Finally, we investigated patient samples based on pCD relapse to identify potential flare predictors ([199]Figure S11). We found that Anaerostipes was significantly depleted in patients experiencing relapse, partially aligning with previous reports of its reduction in IBD patients.[200]^68^,[201]^69^,[202]^70 Anaerostipes is a common gut commensal known for producing butyrate, which is crucial for maintaining the epithelial cell barrier and immune homeostasis.[203]^71^,[204]^72 While recent research has indicated that one of its representative species, Anaerostipes hadrus, may have detrimental effects under certain conditions,[205]^73 our observations of Anaerostipes depletion in relapsed patients suggest that it could serve as a robust biomarker for disease recurrence. Nevertheless, the findings related to Anaerostipes are based on comparisons involving only five patients and thus require further validation in larger cohorts to confirm their robustness and generalizability. This study has certain limitations, such as differences in treatment regimens among patients; specifically, four of the 27 patients did not receive IFX and some of patients did not receive EEN. However, majority of patients ultimately received IFX, and post-treatment samples were collected approximately one year after therapy initiation. Also, note that the potential effects of EEN are assumed to be limited, as the samples were collected approximately 10 months after the completion of the EEN regimen, during which patients had resumed a regular diet. Indeed, when the effects of EEN and IFX were evaluated in relation to changes in both the microbiome and clinical indices ([206]Table S4; [207]Figure S6), their impact appeared to be limited. This suggests that the treatment regimen may not be a primary driver of the observed microbiome alterations and thus was not considered a major factor in our analysis. However, since this study did not incorporate detailed dietary intake assessments, given the dynamic and individualized dietary regimens often followed by pCD patients, future studies incorporating standardized and longitudinal dietary tracking will be essential to better understand the interplay between diet and microbiome dynamics in this population. Another notable limitation of this study is the relatively small sample size, which may have limited the statistical power to detect some associations. Nevertheless, acquiring both stool and paired pre- and post-treatment biopsy samples from pediatric patients poses significant ethical and logical challenges, particularly when targeting multiple GI tract sites. Despite these constraints, our dataset offers a rare and valuable insight into site-specific microbiome and mycobiome dynamics in pCD. To move closer to precision medicine in pCD, future research must be grounded in large-scale, multi-center cohort studies incorporating comprehensive longitudinal sampling of blood, tissue, and stool at critical clinical time points (i.e., before, during, and after treatment).[208]^74 Such efforts are essential not only to validate our current findings, but also to enable the development of reliable biomarkers that reflect disease activity, therapeutic response, and relapse risk. Importantly, the integration of clinical metadata with multi-omics datasets (e.g., microbiome, mycobiome, host transcriptome, and metabolome) collected over time will allow for the systematic tracking of patient trajectories. This framework will facilitate the identification of robust, clinically actionable biomarkers and ultimately support the implementation of personalized treatment strategies in pediatric IBD. Limitations of the study This study offers valuable insight into the spatial and temporal dynamic of the pediatric Crohn’s disease (pCD) microbiome and mycobiome, as well as their clinical correlates before and after treatment. Notably, the simultaneous profiling of both mucosal and luminal bacterial and fungal communities across multiple GI sites in post-treatment samples represents a meaningful contribution to the current understanding of pCD, given the paucity of such spatially resolved analyses in pediatric cohort. Nevertheless, several limitations should be acknowledged when interpreting the findings. First, while the paired integration of stool and biopsy samples provides a valuable longitudinal perspective, the relatively modest cohort size, especially for relapse stratification, may constrain statistical power and limit the generalizability of some associations. This limitation is primarily attributable to the logistical and ethical challenges of repeated mucosal sampling in pediatric populations. Second, although the majority of patients received a standardized treatment protocol, variation in therapeutic exposure (e.g., infliximab, EEN) was not formally stratified in downstream analyses due to sample size constraints. However, considering that post-treatment samples were collected approximately one year after therapy initiation, well beyond the immediate treatment window, and that most patients had resumed a regular diet following EEN, our comparative diversity and compositional analyses suggest that these treatment variables, including infliximab and EEN, did not exert significant effects on the overall microbial landscape. Third, dietary intake, an important modulator of gut microbiota, was not longitudinally or systematically captured, restricting our ability to control for nutritional confounders. Fourth, the absence of healthy control biopsy samples for fungal (ITS2) sequencing limits direct comparisons of site-specific mycobiome alterations. Lastly, the reliance on amplicon-based sequencing (16S rRNA and ITS2) imposes constraints on taxonomic resolution and precludes functional inference beyond predictive modeling, though we sought to mitigate this through careful cross-validation using multiple statistical methods. Resource availability Lead contact Requests for further information and resources should be directed to and will be fulfilled by the lead contact, Dong-Yup Lee (dongyuplee@skku.edu). Material availability This study did not generate new unique reagents. Data and code availability * • Data: The raw 16S rRNA and ITS2 amplicon sequencing data (FASTQ files) generated in this study have been deposited in the NCBI SRA (BioProjects: [209]PRJNA1156939, [210]PRJNA1156940, and [211]PRJNA1156941). * • Code: R scripts used to analyze the data and generate the manuscript figures are available upon request. * • All other items: Any additional information required to reanalyze the data reported in this paper is available from the [212]lead contact upon request. Acknowledgments