Abstract
The infant gut microbiome is essential for long-term health and is linked to atopic dermatitis (AD), although the underlying mechanisms are not fully understood. This study investigated gut microbiome-host interactions in 31 infants with AD and 29 healthy controls using multi-omics approaches, including metagenomic, host transcriptomic, and metabolomic analyses. Microbial diversity was significantly altered in AD, with Bifidobacterium longum and Clostridium innocuum associated with these changes. At the strain-level, only B. longum differed significantly between groups, with pangenome analyses identifying genetic variations potentially affecting amino acid and lipid metabolites. Notably, B. longum subclade I, which was more prevalent in healthy controls, correlated with host transcriptomic pathways involved in phosphatidylinositol 3-kinase-AKT signaling and neuroactive ligand-receptor pathways, as well as specific metabolites, including tetrahydrocortisol and ornithine. These findings highlight the role of B. longum strain-level variation in infants, offering new insights into microbiome-host interactions related to AD.
Similar content being viewed by others
Introduction
Atopic dermatitis (AD), a prevalent, chronic inflammatory skin disorder, affects approximately 20% of children globally, often continuing into adolescence and adulthood1,2. Its development is characterized by a complex and multifactorial pathogenesis involving genetic predisposition, innate and adaptive immune response, compromised skin barrier function, and environmental exposures3. The hygiene hypothesis specifically suggests that reduced microbial exposure in early life contributes to the increasing incidence of AD. Recent studies have highlighted the significance of the gut microbiome and its metabolic activities in relation to human health, strongly associated with several allergic diseases, including asthma and AD4. Furthermore, the Th1/Th2 response of the adaptive immune system is influenced by the gut microbiome, with the Th2 response being particularly dominant in AD and associated with increased immunoglobulin E (IgE) synthesis5.
Given that the physiological systems of an infant are not fully developed until after birth, the first few years of life are critical. During this period, the gut microbiome plays a crucial role in the maturation of both the central nervous system and immune system, alongside the regulation of the stress response, all of which can have lifelong implications6. Recent studies have identified a significant association between the gut microbiome and neurobehavioral outcomes, such as anxiety, depression7, and autism spectrum disorders8, further increasing interest in gut-brain communication9. Furthermore, clinical reports often associate these neurobehavioral findings with skin inflammation, suggesting that the gut-brain-skin axis is a vital communication pathway influenced by neurotransmitters modulated by the gut microbiota10,11,12. The probiotic Bifidobacterium longum has been shown to effectively modulate central nervous system function in both animal models13 and human studies14 while alleviating AD through interactions along the gut-skin axis associated with tryptophan metabolism15. Bifidobacterium is one of the predominant microorganisms in the infant gut microbiome, making the establishment of a symbiotic relationship with healthy Bifidobacterium crucial for future neurodevelopment and skin health.
Recent analyses of the microbiome at the strain-level have gained prominence in understanding human health, including aspects of infection, immunity, nutrition, and disease16,17. The metabolic capabilities of microbial species vary significantly among strains based on their genetic composition16. Consequently, understanding microbial communities with strain-level variations and the associated metabolome is essential for elucidating interactions with the host immune system. Therefore, this study aims to investigate the host–microbial crosstalk associated with AD in infants through comprehensive strain-level resolution metagenomics, metabolomics, and host transcriptomics at 6 months of age. Exploiting these multi-omics aspects could offer insights into the complex relationship between host and microbiome associated with AD.
Results
Effect of the infant gut microbial diversity on atopic dermatitis
A metagenomic analysis was conducted on 60 fecal samples, including 31 from individuals diagnosed with AD and 29 from healthy controls. This analysis incorporated clinical variables, including feeding type, mode of delivery, and family history (Supplementary Table 1 and 2). No significant differences were observed between the groups concerning clinical variables, except for the immunological markers related to the number of samples (p > 0.05, Supplementary Table 1) and the number of sequence reads (p = 0.941, mean = 7,158,967 reads, Supplementary Table 3). These findings support the continuation of further analyses. The predominant microbial species identified by mean relative abundance across all samples was B. longum (37.6%), followed by B. bifidum (15.0%), Escherichia coli (10.6%), and Veillonella parvula (5.9%) (Fig. 1A). However, B. longum was the only species detected in all samples (relative abundance > 0.01%, Fig. 1A).
A Heatmap showing the distribution of the top 25 enriched species between AD patients and healthy controls. B PCoA plot illustrating β-diversity (unweighted UniFrac distance) to compare microbial communities between patients with AD and healthy controls. C Comparison of α-diversity between AD patients and healthy controls. D Pearson correlation between microbial diversity and microbial species (FDR < 0.05). E Structural equation model to assess the mediating effect of microbial diversity on the causal role of B. longum in the AD development (*p < 0.05, **p < 0.01, ***p < 0.001). AD atopic dermatitis, PCoA Principal Coordinate Analysis, PC1 principal component 1.
Although no significant differences in the relative abundance of microbial species were observed between the AD and healthy groups (Wilcoxon rank sum; p > 0.1), significant differences in microbial diversity were observed. The microbial species composition (β-diversity), assessed using unweighted UniFrac distance, revealed substantial differences between groups (Fig. 1B, PERMANOVA; p = 0.023), with principal component 1 (PC1) explaining 22% of the variance in taxonomic composition. Furthermore, α-diversity richness (the number of observed species) was significantly higher in healthy controls compared to infants with AD (Fig. 1C, Wilcoxon rank sum; p < 0.05).
The significant association between microbial diversity and AD prompted an investigation into the specific microbial species influencing this diversity. The analysis revealed a strong correlation between microbial richness and PC1 values (Pearson correlation; p < 10−10), with several microbial species exhibiting correlations with both α- and β-diversity indices, particularly with respect to richness and PC1 values. The relative abundance of B. longum, C. innocuum, and Erysipelatoclostridium ramosum exhibited a significant correlation with both diversity indices (Fig. 1D, false discovery rates (FDR) < 0.05). In contrast, B. breve, Flavonifractor plautii, and Enterococcus faecalis were correlated solely with the PC1 value (Fig. 1D, FDR < 0.05).
Subsequently, formal mediation analyses were conducted to investigate the effects of microbial species on AD through microbial diversity. The mediation analyses indicated that B. longum, C. innocuum, and E. ramosum were not directly related to AD. However, these taxa were associated with variations in microbial diversity, which in turn was associated with AD (Fig. 1E). Among these, B. longum and C. innocuum exhibited full mediation by microbial diversity, demonstrating significant indirect effects on AD. Notably, C. innocuum was fully mediated by both diversity indices (Fig. 1E). While the relative abundance of B. longum and C. innocuum did not significantly differ between the AD group and healthy controls, both species showed a strong association with microbial diversity, suggesting a potential indirect role in AD (Fig. 1E, p < 0.05).
Distinct subclades of Bifidobacterium longum subspecies infantis are enriched in different skin phenotypes
To investigate the colonization of infants with AD and healthy controls through different bacterial strains, a strain-level analysis was conducted on three predominant species: B. longum, E. coli, and V. parvula. These species were adequately sequenced across all groups (Fig. 1A). Among these, only B. longum exhibited strain-level stratification based on skin phenotype, with subclades I and II associated with healthy controls and infants with AD, respectively (Fisher’s exact test, p < 0.05; Fig. 2A and Supplementary Fig. 1). Both subclades were classified as B. longum subsp. infantis, while other strains showed no specific association with disease status.
A Phylogenetic analysis based on StrainPhlan3, derived from marker genes of SNVs in B. longum, reveals two distinct subclades corresponding to AD and healthy controls, along with the pangenome distribution of KO via PanPhlan3. Feeding types are indicated with superscripts: b for breastfeeding, f for formula feeding, and no superscript for mixed formula. B Comparison of ANI values among B. longum MAGs, showing cluster patterns within and between subclades. C PCoA plot based on pangenome clusters generated using Panaroo, illustrating differences in genetic profiles between subclades and AD. AD atopic dermatitis, KO KEGG orthologs, ANI Average Nucleotide Identity.
To validate strain differentiation, B. longum metagenome-assembled genomes (MAGs) were reconstructed for each individual, yielding high-quality and medium-quality MAGs (mean completeness: 93.3%, contamination: 2.6%, Supplementary Table 4, https://doi.org/10.6084/m9.figshare.27367887) in 55 out of 60 individuals (91.7%). The reconstructed MAGs were classified into B. longum subsp. infantis subclades I and II. Phylogenetic analyses and average nucleotide identity (ANI) comparisons differentiated these subclades within B. longum subsp. infantis (Fig. 2B). Subsequent pangenome analysis revealed distinct gene distributions across the subclades, suggesting differential functional potential.
Overall, 258 KEGG orthologs (KOs) were evaluated using PanPhlan analysis to identify genes enriched in each subclade of B. longum subsp. infantis. Twenty-four genes exhibited significant differences between subclades I and II (FDR < 0.05; Fisher’s exact test), indicating distinct functional profiles compared to other B. longum subsp. infantis strains outside these subclades (Fig. 2A). Genes associated with bacterial defense mechanisms showed variation between the subclades, with subclade I possessing the CRISPR-associated protein Cas1 (K15342), Cas3 (K07012), Cas5d (K19119) and subclade II containing elements of the Restriction-Modification (RM) system (K03427) along with a toxin gene (K06218). Furthermore, three genes were identified as potential contributors to amino acid metabolism. The gene specific to subclade I encoded an asparagine synthase (glutamine-hydrolyzing enzyme; K01953), which catalyzes the conversion of aspartate to asparagine while converting glutamine to glutamate. In contrast, the subclade II pangenomes contained a proline/betaine transporter (K03762), a glutamate transport system (K10008), a branched-chain amino acid (BCAA) transporter (K01997, K01998), an amidohydrolase (metallo peptidase: K01436), and dipeptidyl-peptidase IV (K01278).
MAG-based pangenome analysis was conducted to capture strain-level variations more precisely than reference-based methods. Through this approach, 777 gene clusters with significant differentiation between subclades were identified (FDR < 0.1), with 16 of these gene clusters showing significant associations with AD status and the healthy control group (Supplementary Table 5). Principal Coordinate Analysis (PCoA) based on the presence or absence of these 16 genes demonstrated significant differences associated with AD status and B. longum subclades (Fig. 2C, PERMANOVA; B. longum type: p = 0.001, AD: p = 0.002). Most gene variations observed were related to protein (dap4, dipeptidyl-peptidase IV; p = 0.0001), amino acid (proC, Pyrroline-5-carboxylate reductase; p = 0.0006), lipid (group_1991, short-chain dehydrogenase; p = 0.0004), and lactose metabolism (purR, lactose operon repressor; p = 0.0005), highlighting pathways potentially influencing host metabolic processes. Additionally, variations in the serine/threonine-protein kinase gene (pknB; p < 0.0002) and forkhead-associated ___domain-containing genes (group_104, group_1440, group_220; p < 0.0005) may differentially influence bacterial signaling pathways across B. longum subclades by modulating protein interactions. In the reference-based PanPhlan analysis, dap4 was prevalent only in subclade II, while the MAG-based pangenome analysis confirmed that different gene clusters of dap4 were present between each clade, suggesting potential differences in peptidase substrates between the two clades.
Multi-omics analysis reveals potential crosstalk between Bifidobacterium longum subclades and host interactions
The colonocyte transcriptome associated with B. longum was further examined based on its subclades. To minimize confounding effects, feeding type, mode of delivery, and family history were adjusted as fixed effects before calculating the correlation between B. longum relative abundance and host transcripts. Overall, 71 and 53 transcripts were significantly correlated with clade I and II, respectively (Fig. 3, p < 0.01, |r | > 0.5). Enrichment analysis was performed on these transcripts to identify the pathways involved (FDR < 0.05). For subclade I, most transcripts were associated with multiple pathways; however, a distinct positive correlation was observed with the neuroactive ligand-receptor interaction pathway (GABRR2, GPR156, HTR1B, GRM4, ADRA2B). Additionally, subclade I exhibited a negative correlation with monoamine oxidase B (MAOB), an enzyme that catalyzes the breakdown of monoamine neurotransmitters into inactive metabolites18, suggesting a potential link to the gut-brain axis. Subclade I was correlated with transcripts involved in several cellular processes, including autophagy, mammalian target of rapamycin (mTOR) signaling, and T cell receptor signaling pathway (AKT2, TNF, MAPK9). In contrast, subclade II was associated with transcripts linked to the longevity-regulating pathway, while other transcripts exhibited no significant enriched pathways (FDR > 0.05).
Correlation network showing the relationships between B. longum subsp. infantis subclades, host colonocyte transcripts, and gut metabolites. Octagonal, circular, rectangular, and triangular nodes represent B. longum subclades, colonocyte transcripts, pathway names, and metabolites, respectively. Nodes connected to subclades are highlighted with a thick outline. Edges indicate significant Spearman correlations (p < 0.01 and |r | > 0.5; positive: blue and negative: red) between the residuals of species, transcripts, and metabolites, adjusted by a generalized linear regression model using feeding type, delivery mode and family history as fixed effects.
Further untargeted metabolomic analyses were performed on fecal samples classified by the B. longum subclade group for additional interpretation. However, the nature of the cohort limited the availability of identical samples for further study (subclade I: 8, subclade II: 17). This analysis identified five and 15 metabolites differentially correlated with subclades I and II, respectively (Fig. 3, Supplementary Table 6, p < 0.01, |r | > 0.5), with no overlap between the metabolites of the two subclades. Metabolites were correlated with the human transcriptome; however, no significant correlations were identified that warranted integration of all omics results. The only metabolite found to correlate with both the human transcriptome and subclade I was m/z 132.02 (1-benzothiophene).
Some metabolite associations appeared to be influenced by pangenomic differences between subclades. Most metabolites associated with the subclades were related to lipid metabolism, with correlation directions differing by subclade. Subclade II exhibited positive correlations with nine lipid-related metabolites (Fig. 3, Supplementary Table 6), while subclade I demonstrated negative correlations with m/z 294.22 (13-L-hydroperoxylinoleic acid, tetrahydrocortisol, 13-oxoODE, 9(S)-HPODE, and stearidonic acid). These findings suggest that the presence of the short-chain dehydrogenase gene (group_1991) in subclade II is positively associated with various lipid-related metabolites, while subclade I exhibits fewer lipid metabolites and demonstrates negative correlations. Furthermore, subclade I showed a positive correlation with the only metabolite related to amino acid metabolism, m/z 132.10 (ornithine, 2,4-diaminopentanoate). Ornithine serves as a precursor to glutamate and proline, with this association potentially attributed to genetic differences in K10008, K01953, K03762 and proC between the subclades.
Discussion
This study reevaluates the role of Bifidobacterium in the gut microbiome of 6-month-old infants, emphasizing its association with AD. While previous research indicated that variability in Bifidobacterium is influenced by feeding type19, this analysis has been expanded to include species-, and strain-level differentiation. Specific strains of B. longum subsp. infantis exhibit distinct phylogenetic and genetic patterns of colonization in individuals with AD compared to healthy individuals, independent of feeding type (Fig. 2A).
Previous studies investigating the relationship between the infant microbiome and AD have primarily relied on 16S rRNA sequencing. However, this approach has some limitations, including reduced taxonomic resolution20, the potential for diversity inflation21, and bias issues22 during amplification. Furthermore, host interactions have been inferred primarily from gut microbiome profiles and metabolic pathway analyses19, without incorporating multi-omics data to directly assess crosstalk, thereby limiting the understanding of AD development. To address these limitations, strain-level metagenomic analysis was used to compare healthy infants and those with AD, integrating findings with metabolomics and host colonocyte transcriptomics to clarify host associations. The findings suggest that B. longum is a major species associated with microbial diversity, with the enrichment of specific B. longum subsp. infantis strains closely related to host interactions. These changes appear most pronounced within the immune and nervous systems of the gut, suggesting they may be mediated through shifts in strain-specific metabolite profiles. Given that disruptions in the infant gut microbiome can significantly affect early biological, immune, and psychological development, these findings hold particular relevance for understanding allergic diseases.
Bifidobacterium is widely recognized as a beneficial probiotic; however, its association with AD has historically been overlooked19,23. Recent reports24,25 indicate an overrepresentation of Bifidobacterium in the gut microbiome of children, which may correlate with allergies. Considering that even beneficial gut microbes can occasionally induce disease-associated imbalances26, the findings suggest that B. longum may be indirectly associated with AD through its modulation of microbial diversity. In this study, AD-associated diversification, specifically at the strain-level within B. longum, was identified in this study. However, additional mediation analysis revealed that B. longum strain type could directly explain AD directly (total & direct effect, p < 0.001) rather than mediating through microbial diversity (indirect [diversity] effect, p = 0.377). Subsequent pangenome and multi-omics analysis further suggests that varying colonization strategies among subclades may influence the different host response directly associated with AD.
The two B. longum subclades exhibited distinct bacterial defense mechanisms, encompassing CRISPR, RM, and toxin-antitoxin (TA) systems27. These variations suggest potential differences in phage susceptibility and host range27, which could be pivotal for niche-specific adaptations within the infant gut microbiome. Subclade II contains a toxin gene (K06218) associated with the type II TA system. TA systems regulate host cell death, the formation of persister states28, and responses to environmental changes28,29, which could represent another genetic determinant influencing colonization strategies.
Infants, with their still-developing hypothalamic-pituitary-adrenal axis and gut microbiome, demonstrate increased sensitivity to stressors. This vulnerability can disrupt the composition of their gut microbiome, potentially leading to conditions such as leaky gut syndrome30. The examination of healthy infants enriched with B. longum subclade I revealed correlations with the phosphatidylinositol 3-kinase (PI3K)-Akt signaling pathway in the host transcriptome, further complemented by associations with the MAPK pathway. These signaling cascades play a critical role in regulating intestinal epithelial cell proliferation, differentiation, and survival31. Furthermore, they stimulate the mTOR, modulating various cellular metabolic processes. TNF within intestinal epithelial cells can activate both survival and apoptotic pathways32. Additionally, TNF can induce PI3K-AKT expression in response to survival signaling, with cell fate dependent on the specific cellular context32. These findings are consistent with those of previous studies that the PI3K-Akt pathway is more pronounced in the microbiome metabolic pathways of healthy infants than in those with AD19. These findings highlight that the expression of transcripts involved in this pathway is associated with the colonocyte of the host, suggesting that specific strains of B. longum may modulate the vitality and activity of intestinal epithelial cells. The increased relative abundance of subtype I, alongside the decreased levels of tetrahydrocortisol and increased transcripts for TNF, AKT2, and MAPK9 in host colonocytes, appears to positively influence intestinal immune maturation.
BCAAs are essential for neonatal nutrition, supplying a significant portion of the vital protein requirements33,34. They are involved in various metabolic functions, such as glucose metabolism, the promotion of the innate immune system, and the suppression of harmful gut microbes35,36. Additionally, BCAAs, along with glutamate—an essential excitatory neurotransmitter37— directly and indirectly affect brain function by stimulating the synthesis of aromatic amino-acid-based neurotransmitters33,38. Hormones in the gut significantly influence the gut-brain axis via the vagus nerve39 and can be influenced by genes such as bacterial dap4 and host MAOB. In accordance with evidence indicating that different Bifidobacterium strains can differentially regulate the gut-brain axis40,41,42, this study suggests that genetic variations among B. longum strains, alongside corresponding differences in host transcript expression, may lead to alterations in gut hormones and neurotransmitter levels, ultimately resulting in distinct regulatory effects on the gut-brain axis. B. longum infantis subclade I exhibits genetic differences that may influence metabolite production, specifically affecting the enteric nervous system. This subclade has been positively correlated with ornithine and associated with various neurotransmitter receptor transcripts in the host. This finding suggests that distinct commensal bacteria produce hormones30 that activate the enteric nervous system and stabilize the gut-brain axis in response to diverse signals.
Psychological factors are closely linked to allergic conditions such as asthma and AD, suggesting potential bidirectional relationships3,43. Different Bifidobacterium strains can affect various aspects of microbiome development44, which may extend to gut-brain axis regulation. These findings highlight the need for further research to clarify if and how psychological factors during early childhood contribute to AD development.
In this study, the small sample sizes from subclades in each group restricted the performance of relevant multi-omics analyses. Furthermore, due to the complex composition of feces, analyzing fecal metabolites may yield inconsistent interpretations. Consequently, although it was challenging to identify correlations across all three omics layers (metagenome, human transcriptome, and metabolome), the genomic and metabolomic analysis of subclade I—characteristic of a healthy population—followed by interpretation of the human transcriptome, demonstrates the potential for meaningful biological insights. This approach underscores the value of a comprehensive, multi-omics strategy for uncovering novel insights.
In essence, we identified significant differences in microbial diversity between healthy infants and those with AD, with dominant B. longum strains playing key roles. Furthermore, distinct B. longum strains clustered into two separate subclades, corresponding to either infants with AD or healthy infants, though this requires validation in larger, independent cohorts. These subclades exhibited unique genetic profiles, metabolite associations, and correlations with host transcriptomes. Our findings underscore the significant influence of specific B. longum strains on gut immune development and the stability of the enteric nervous system. Additionally, early colonization by these strains may be crucial in shaping the skin health of the host.
Methods
Participants and clinical assessments
The study included patients from the Cohort for Childhood Origin of Asthma and Allergic Diseases, a longitudinal, general population-based birth cohort, and the Childhood Asthma Atopy Center at Asan Medical Center22. The study population consisted of sixty 6-month-old infants, including 29 healthy controls and 31 infants with AD. Baseline characteristics are provided in Supplementary Table 1 and 2. This study was conducted in accordance with the Declaration of Helsinki and approved by the Institutional Review Board (IRB) of Asan Medical Center (IRB No. 2008-0616 and 2015-1031). Written informed consent was obtained from the parents of all participants.
Pediatric allergists diagnosed AD based on the criteria of Hanifin and Rajka45. The severity of AD was evaluated using the Scoring Atopic Dermatitis index. Total and specific serum IgE levels (for egg white and milk)(IU/mL) were measured using the ImmunoCAP-CAP 1000 system (Phadia AB, Uppsala, Sweden). Blood eosinophil percentages were measured at 6 months of age using an automatic blood cell counter (XE-100; Sysmex Co., Kobe, Japan).
Preprocessing of metagenomic data
Metagenomic sequencing was performed on the gut microbiota from 60 fecal samples of 6-month-old infants. DNA was extracted using the RNeasy PowerMicrobiome Kit (Qiagen, Valencia, CA, USA) and then fragmented with NEBNext dsDNA Fragmentase (Cat #0348 L, New England Biolabs, Ipswich, MA, USA). The sequencing library was prepared using ACCEL-NGS 2S PLUS DNA Library Kits (Cat #21096, Swift Biosciences, Ann Arbor, MI, USA) following the instructions of the manufacturer. After verifying the metagenomic library size with a Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA), we prepared the library as previously described by Lee et al. 22. Equimolar concentrations (2 nM) of each library were quantified through quantitative real-time PCR using a TaKaRa PCR Thermal Cycler Dice Real Time System III (TaKaRa Bio, Inc., Shiga, Japan) and the GenNext NGS Library Quantification Kit (Cat #NLQ-101, Toyobo, Osaka, Japan). Sequencing was conducted using the Illumina HiSeq 2500 system, producing 250 bp paired-end reads.
Quality control of shotgun metagenomic raw reads was completed using FaQCs46, using the -q 30 flag in default mode to trim low-quality reads. Following this, host-derived sequences were removed using BWA-MEM47 against the human reference genome GRCh38.p13. We also employed FastUniq48 to correct duplicate paired-end read errors occurring from the Illumina pattern flow cell method. Supplementary Table 3 shows the basic statistics of quality-filtered metagenomic data.
Metagenome profiling
Taxonomic profiling at the species level was conducted using MetaPhlAn349 with default settings (Supplementary Table 7). Microbial diversity was calculated (Supplementary Table 8) using the “calculate_diversity.R” script, available on the MetaPhlAn GitHub repository. To compute unweighted UniFrac distances, we employed the precomputed phylogenetic tree provided by MetaPhlAn3 (mpa_v30_CHOCOPhlAn_201901_species_tree.nwk). The unweighted UniFrac distance was then utilized for subsequent PCoA analysis. A PERMANOVA test comparing AD and healthy control groups was performed using the Vegan package50 in R.
Mediation analysis of the AD phenotype was conducted using the B. longum relative abundance (transformed with an arcsine square root) and microbial diversity, utilizing the lavaan package in R. Additionally, we performed a mediation analysis using B. longum subclades instead of relative abundance to assess the influence of each subclade on AD, as represented by the formula below.
# Direct effect: Disease ~ c*(Blongum_abundance or subclades) + b*diversity (observed species or PC1 value)
# Mediator: diversity (observed species or PC1 value)~ a* (Blongum_abundance or subclades)
# Indirect effect: (a*b)
We utilized StrainPhlAn (default option)49 for strain-level microbiome analysis. Consensus marker sequences were screened from mapped reads to the MetaPhlAn marker database, extracting only the marker sequences of abundant microbial species per sample. B. longum, Escherichia coli, and Veillonella parvula (Supplementary Fig. 1) were the most abundant species across samples. Their marker sequences were used to calculate strain-level phylogenetic comparisons with db_markers provided by MetaPhlAn3. PanPhlAn (default option)49 was also used to identify genetic differences between these strains and independently generated the presence or absence profiles of functional genes for each sample using the reference pangenome database. Finally, genetic differences across lineages were visualized using iTOL51, based on data produced by StrainPhlAn and PanPhlAn.
Pangenomic analysis of B. longum recovered by metagenome-assembled gnomes
Quality-filtered metagenomic reads were assembled into contigs using MEGAHIT52 (v1.2.9; with the meta-sensitive option). Genome binning was subsequently performed with MetaBAT253 (v.2.12.1; default settings), CONCOCT54 (v.1.0.0; default settings), and VAMB55 (v.3.0.3.2; default settings). Further refinement of each genome bin was conducted with ACR56 (v.0.2; default settings) to improve bin quality, followed by a dereplication step using dRep57 (v.3.2.0; options -comp 50 -con 10) to select the highest-quality MAGs from each sample. Taxonomic classification was carried out using the Genome Taxonomy Database Toolkit58 (GTDB-Tk, v.2.0.0; default settings), from which B. longum was selected for subsequent pangenome analysis. Pangenome clusters were identified using Panaroo59 (v1.4.0; with the -clean-mode strict option), and each MAG was compared using fastANI60 (v1.33; default settings) to calculate ANI values.
For additional strain-level pangenome differences, we first filtered the presence of gene clusters occurring in over 25% of all samples. Fisher’s exact test was then conducted to compare the AD and healthy control groups, as well as the B. longum subsp. clades. After applying the subsequent FDR test, we selected significant genes that discriminate between AD and B. longum subsp. clades. The number of gene clusters present in each sample was used to calculate a distance matrix using Bray-Curtis distance. PCoA analysis and PERMANOVA tests were then performed using the Vegan package50 in R.
Human transcriptome processing
Exfoliated colonocytes were isolated from 60 fecal samples and stored at -70°C until analysis, using the Percoll-density gradient centrifugation method, as described in previous studies61,62,63. Briefly, thawed fecal samples (0.5 g) were vortexed with 10 mL of phosphate-buffered saline (PBS) and filtered through a 40 μm cell strainer (SPL, Seoul, Korea). The filtrate was carefully layered onto Histopaque-1077 (Sigma Aldrich, St. Louis, MO) and centrifuged at 400 × g for 30 min at room temperature. The resulting pellets were then washed twice with 10 mL PBS. Total mRNA was isolated from these colonocytes using an RNeasy Mini Kit (Qiagen, Hilden, Germany). The synthesized cDNA was hybridized onto a GeneChip® Human Gene 2.0 ST Array (Affymetrix) per the protocol of the manufacturer, using the GeneChip WT Pico Reagent kit (Affymetrix, Santa Clara, CA). After array scanning with the GCS3000 Scanner (Affymetrix), raw data were normalized using the Robust Multichip Analysis algorithm within Affymetrix Power Tools. As the experiment involved two datasets, batch effect correction was applied using the ComBat package64.
Global metabolome profiling
Thirty-four fecal samples were collected and immediately stored at -80°C until metabolite analysis. Fecal metabolites were extracted using a standard liquid-liquid separation technique65,66. Briefly, a 2:1 mixture of chloroform and methanol was added to the frozen feces, followed by centrifugation for 15 min. Polar metabolites were collected from the upper aqueous phase, while nonpolar, lipid-containing metabolites were obtained from the lower organic phase. Liquid chromatography-mass spectrometry (LC-MS) analysis was performed using an Ultimate 3000 (Dionex) and an LTQ Orbitrap XL (Thermo Fisher). A reverse-phase column (Pursuit 5; 150 × 2.0 mm) was employed for separating nonpolar metabolites, whereas a HILIC column (HILIC Plus; 100 × 2.1 mm) was used for polar metabolites. LC-MS analysis was conducted for each sample solution in positive and negative ion modes. Metabolite features, including mass-to-charge ratios (m/z) and retention times, were extracted using Compound Discoverer 2.0.
Correlations between Bifidobacterium longum subclades and multi-omics data
We constructed a generalized linear model using MaAsLin267 to determine the significant associations between multi-omics and B. longum subclades. In this model, each omics dataset—metagenomic species (60 samples), human transcriptome (60 samples), and metabolome (34 samples) —was adjusted for factors including the feeding type, mode of delivery, family history, and AD phenotype of the infant. Subsequent analyses were performed according to the B. longum subclade (-I: 12 samples and -II: 29 samples), although for the metabolome, only 8 and 17 samples were available for subclade-I and -II, respectively. Pearson residuals from this model were used to evaluate correlations between the relative abundance of each B. longum subclade and the human transcriptome and metabolome. We utilized the Spearman method to determine correlations, selecting features with a p-value of < 0.01, an absolute r value of > 0.5 for host transcripts and fecal metabolites. Gene enrichment analysis was then performed using Enrichr68 to identify pathways significantly associated with the host transcriptome, considering KEGG pathways with an adjusted p-value of < 0.05 as significant. Finally, Cytoscape was used to visualize these identified correlations and significant pathways69.
Data availability
The raw metagenome sequencing data generated during the current study are available in the NCBI SRA under BioProject accession number PRJNA979436 and in the ENA SRA under accession number PRJEB45443 (Supplementary Table 2). B. longum MAGs have been deposited in Figshare under the following: https://doi.org/10.6084/m9.figshare.27367887.
References
Langan S., Irvine A., Weidinger S. Atopic dermatitis LANCET. 396, 758–758 (2020).
Kim, J. P., Chao, L. X., Simpson, E. L. & Silverberg, J. I. Persistence of atopic dermatitis (AD): a systematic review and meta-analysis. J. Am. Acad. Dermatol. 75, 681–687.e611 (2016).
Lee, S.-Y., Lee, E., Park, Y. M. & Hong, S.-J. Microbiome in the gut-skin axis in atopic dermatitis. Allergy, asthma Immunol. Res. 10, 354–362 (2018).
Huang, Y. J. et al. The microbiome in allergic disease: Current understanding and future opportunities-2017 PRACTALL document of the American Academy of Allergy, Asthma & Immunology and the European Academy of Allergy and Clinical Immunology. J. Allergy Clin. Immunol. 139, 1099–1110 (2017).
Ilves, L. et al. Metabolomic analysis of skin biopsies from patients with atopic dermatitis reveals hallmarks of inflammation, disrupted barrier function and oxidative stress. Acta Derm.-Venereologica 101, adv00407–adv00407 (2021).
Derrien, M., Alvarez, A. S. & de Vos, W. M. The Gut Microbiota in the First Decade of Life. Trends Microbiol 27, 997–1010 (2019).
Kelly, J. R., O’Keane, V., Cryan, J. F., Clarke, G. & Dinan, T. G. Mood and microbes: gut to brain communication in depression. Gastroenterol. Clin. 48, 389–405 (2019).
Liu, F. et al. Altered composition and function of intestinal microbiota in autism spectrum disorders: a systematic review. Transl. psychiatry 9, 43 (2019).
Zhu, S. et al. The progress of gut microbiome research related to brain disorders. J. neuroinflammation 17, 1–20 (2020).
Bowe, W. P. & Logan, A. C. Acne vulgaris, probiotics and the gut-brain-skin axis-back to the future?. Gut Pathog. 3, 1–11 (2011).
Wang, X. et al. Dysregulation of the gut-brain-skin axis and key overlapping inflammatory and immune mechanisms of psoriasis and depression. Biomedicine Pharmacother. 137, 111065 (2021).
Chang, H. Y. et al. Prenatal maternal distress affects atopic dermatitis in offspring mediated by oxidative stress. J. Allergy Clin. Immun. 138, 468 (2016).
Luk, B. et al. Postnatal colonization with human” infant-type” Bifidobacterium species alters behavior of adult gnotobiotic mice. PLoS One 13,e0196510 (2018).
Wang, H., Braun, C., Murphy, E. F. & Enck, P. Bifidobacterium longum 1714™ strain modulates brain activity of healthy volunteers during social stress. Am. J. Gastroenterol. 114, 1152 (2019).
Fang, Z. et al. Bifidobacterium longum mediated tryptophan metabolism to improve atopic dermatitis via the gut-skin axis. Gut Microbes 14, 2044723 (2022).
Park, S. Y., et al. Strain-level fitness in the gut microbiome is an emergent property of glycans and a single metabolite. Cell 185, 513 (2022).
Chen Y. W., et al. Gut metagenomes of type 2 diabetic patients have characteristic single-nucleotide polymorphism distribution in Bacteroides coprocola. Microbiome. 5. (2017)
Zhu, J. X., ed. Dopamine in the Gut. Springer, 2021.
Lee, M. J. et al. Perturbations of gut microbiome genes in infants with atopic dermatitis according to feeding type. J. Allergy Clin. Immunol. 141, 1310–1319 (2018).
Hildebrand, F. Ultra-resolution Metagenomics: When Enough Is Not Enough. mSystems 6, e0088121 (2021).
Schloss, P. D. Amplicon Sequence Variants Artificially Split Bacterial Genomes into Separate Clusters. Msphere 6, e0019121 (2021).
Lee, M. J. et al. Disordered development of gut microbiome interferes with the establishment of the gut ecosystem during early childhood with atopic dermatitis. Gut Microbes 14, 2068366 (2022).
Ta, L. D. H. et al. A compromised developmental trajectory of the infant gut microbiome and metabolome in atopic eczema. Gut microbes 12, 1801964 (2020).
Simonyté Sjödin, K. et al. Temporal and long-term gut microbiota variation in allergic disease: a prospective study from infancy to school age. Allergy 74, 176–185 (2019).
Štšepetova, J. et al. Molecularly assessed shifts of Bifidobacterium ssp. and less diverse microbial communities are characteristic of 5-year-old allergic children. FEMS Immunol. Med. Microbiol. 51, 260–269 (2007).
Wang, W. et al. Increased proportions of Bifidobacterium and the Lactobacillus group and loss of butyrate-producing bacteria in inflammatory bowel disease. J. Clin. Microbiol 52, 398–406 (2014).
Bernheim, A. & Sorek, R. The pan-immune system of bacteria: antiviral defence as a community resource. Nat. Rev. Microbiol. 18, 113–119 (2020).
Wen, Y., Behiels, E. & Devreese, B. Toxin–Antitoxin systems: their role in persistence, biofilm formation, and pathogenicity. Pathog. Dis. 70, 240–249 (2014).
Zhang, L.-Y. et al. Toxin-antitoxin systems alter adaptation of Mycobacterium smegmatis to environmental stress. Microbiol. Spectr. 10, e02815–e02822 (2022).
Sarkar, A. et al. Psychobiotics and the Manipulation of Bacteria-Gut-Brain Signals. Trends Neurosci. 39, 763–781 (2016).
Mandal, C. C., Ghosh Choudhury, G. & Ghosh-Choudhury, N. Phosphatidylinositol 3 kinase/Akt signal relay cooperates with smad in bone morphogenetic protein-2-induced colony stimulating factor-1 (CSF-1) expression and osteoclast differentiation. Endocrinology 150, 4989–4998 (2009).
Yan, F., John, S. K. & Polk, D. B. Kinase suppressor of Ras determines survival of intestinal epithelial cells exposed to tumor necrosis factor. Cancer Res. 61, 8668–8675 (2001).
Amari, S., Shahrook, S., Namba, F., Ota, E. & Mori, R. Branched-chain amino acid supplementation for improving growth and development in term and preterm neonates. Cochrane Db. Syst. Rev. 10, CD012273 (2020).
de Groof, F. et al. Branched-chain amino acid requirements for enterally fed term neonates in the first month of life. Am. J. Clin. Nutr. 99, 62–70 (2014).
Layman, D. K. & Walker, D. A. Potential importance of leucine in treatment of obesity and the metabolic syndrome. J. Nutr. 136, 319s–323s (2006).
Alam, N. H. et al. L-isoleucine-supplemented Oral Rehydration Solution in the Treatment of Acute Diarrhoea in Children: A Randomized Controlled Trial. J. Health Popul Nutr. 29, 183-190 (2011).
Filpa, V. et al. Role of glutamatergic neurotransmission in the enteric nervous system and brain-gut axis in health and disease. Neuropharmacology 111, 14–33 (2016).
Fernstrom, J. D. Branched-chain amino acids and brain function. J. Nutr. 135, 1539S–1546S (2005).
Matsubara, Y., Kiyohara, H., Teratani, T., Mikami, Y. & Kanai, T. Organ and brain crosstalk: The liver-brain axis in gastrointestinal, liver, and pancreatic diseases. Neuropharmacology 205, 108915 (2022).
Tian, P., Wang, G., Zhao, J., Zhang, H. & Chen, W. Bifidobacterium with the role of 5-hydroxytryptophan synthesis regulation alleviates the symptom of depression and related microbiota dysbiosis. J. Nutr. Biochem 66, 43–51 (2019).
Strandwitz, P. Neurotransmitter modulation by the gut microbiota. Brain Res. 1693, 128–133 (2018).
Savignac, H. M., Kiely, B., Dinan, T. G. & Cryan, J. F. Bifidobacteria exert strain-specific effects on stress-related behavior and physiology in BALB/c mice. Neurogastroent Motil. 26, 1615–1627 (2014).
Chida, Y., Hamer, M. & Steptoe, A. A bidirectional relationship between psychosocial factors and atopic disorders: A systematic review and meta-analysis. Psychosom. Med. 70, 102–116 (2008).
Ma, B. et al. Strain-specific alterations in gut microbiome and host immune responses elicited by tolerogenic Bifidobacterium pseudolongum. Sci. Rep. 13, 6628 (2023).
Hanifin, J. M. Diagnostic features of atopic dermatitis. Acta Derm. Venereol. 92, 44–47 (1980).
Lo, C.-C. & Chain, P. S. Rapid evaluation and quality control of next generation sequencing data with FaQCs. BMC Bioinform. 15, 1–8 (2014).
Li, H. Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. arXiv preprint arXiv:13033997. (2013).
Xu, H. et al. FastUniq: a fast de novo duplicates removal tool for paired short reads. PloS one 7,e52249 (2012).
Beghini, F. et al. Integrating taxonomic, functional, and strain-level profiling of diverse microbial communities with bioBakery 3. elife 10, e65088 (2021).
Dixon, P. VEGAN, a package of R functions for community ecology. J. vegetation Sci. 14, 927–930 (2023).
Letunic, I. & Bork, P. Interactive Tree Of Life (iTOL) v5: an online tool for phylogenetic tree display and annotation. Nucleic acids Res. 49, W293–W296 (2021).
Li, D. et al. MEGAHIT v1. 0: a fast and scalable metagenome assembler driven by advanced methodologies and community practices. Methods 102, 3–11 (2016).
Kang, D. D. et al. MetaBAT 2: an adaptive binning algorithm for robust and efficient genome reconstruction from metagenome assemblies. PeerJ. 7, e7359 (2019).
Alneberg, J. et al. Binning metagenomic contigs by coverage and composition. Nat. Methods. 11, 1144–1146 (2014).
Nissen, J. N. et al. Improved metagenome binning and assembly using deep variational autoencoders. Nat. Biotechnol. 39, 555–560 (2021).
Seong, H. J., Kim, J. J. & Sul, W. J. ACR: metagenome-assembled prokaryotic and eukaryotic genome refinement tool. Brief. Bioinform. 24, bbad381 (2023).
Olm, M. R., Brown, C. T., Brooks, B. & Banfield, J. F. dRep: a tool for fast and accurate genomic comparisons that enables improved genome recovery from metagenomes through de-replication. ISME J. 11, 2864–2868 (2017).
Chaumeil, P. A. et al. GTDB-Tk: a toolkit to classify genomes with the Genome Taxonomy Database. Bioinformatics 36, 1925–1927 (2020).
Gerry, T. H. et al. Producing polished prokaryotic pangenomes with the Panaroo pipeline. Genome Biol. 21, 1–21 (2020).
Chirag, J. et al. High throughput ANI analysis of 90K prokaryotic genomes reveals clear species boundaries. Nat. Commun. 9, 5114 (2018).
Albaugh, G. P. et al. Isolation of exfoliated colonic epithelial cells, a novel, non-invasive approach to the study of cellular markers. Int. J. Cancer 52, 347–350 (1992).
Chandel, D. S., Braileanu, G. T., Chen, J.-H. J., Chen, H. H. & Panigrahi, P. Live colonocytes in newborn stool: surrogates for evaluation of gut physiology and disease pathogenesis. Pediatr. Res. 70, 153–158 (2011).
Yamao, T. et al. Abnormal expression of CD44 variants in the exfoliated cells in the feces of patients with colorectal cancer. Gastroenterology 114, 1196–1205 (1998).
Leek, J. T., Johnson, W. E., Parker, H. S., Jaffe, A. E. & Storey, J. D. The sva package for removing batch effects and other unwanted variation in high-throughput experiments. Bioinformatics 28, 882–883 (2012).
Bligh, E. G. & Dyer, W. J. A rapid method of total lipid extraction and purification. Can. J. Biochem Physiol. 37, 911–917 (1959).
Folch, J., Lees, M. & Sloane Stanley, G. H. A simple method for the isolation and purification of total lipides from animal tissues. J. Biol. Chem. 226, 497–509 (1957).
Mallick, H. et al. Multivariable association discovery in population-scale meta-omics studies. PLoS computational Biol. 17,e1009442 (2021).
Chen, E. Y. et al. Enrichr: interactive and collaborative HTML5 gene list enrichment analysis tool. BMC Bioinform. 14, 1–14 (2013).
Shannon, P. et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 13, 2498–2504 (2003).
Acknowledgements
This study was supported by Korea National Institute of Health (KNIH) research project (2024ER211600); the Korea Environment Industry & Technology Institute (KEITI) through the Core Technology Development Project for Environmental Diseases Prevention and Management funded by the Ministry of Environment (MOE) (RS-2022-KE002048); the Bio & Medical Technology Development Program of the National Research Foundation of Korea (NRF), funded by the Korean government (MSIT) (NRF-2017M3A9F3043834); and the Korean Centers for Disease Control and Prevention (2008-E33030-00, 2009-E33033-00, 2011-E33021-00, 2012-E33012-00, 2013-E51003-00, 2014-E51004-00, 2014-E51004-01, 2014-E51004-02, 2017-E67002-00, 2017-E67002-01, 2017-E67002-02, 2020E670200, 2020E670201 and 2020E670202).
Author information
Authors and Affiliations
Contributions
Methodology, H.J.S., J.H.K., and H.J.Y.; Visulization, H.J.S.; Managing Omics Data, Y.M.P., T.Y.K., and S.M.Y.; Writing – Original Draft, H.J.S., Y.M.P., and S.J.H.; Writing – Review and Editing, H.J.S., Y.M.P., J.H.K. and S.J.H.; Cohort Data Management, S.Y.L.; Microbiome Data Management, Y.K.L., B.S.K., D.W.L. and M.H.N.; Funding Acquisition, S.J.H.; Supervision, S.J.H.
Corresponding author
Ethics declarations
Competing interests
The authors declare no competing interests.
Additional information
Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary information
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.
About this article
Cite this article
Seong, H.J., Park, Y.M., Kim, BS. et al. Integrated multi-omics reveals different host crosstalk of atopic dermatitis-enriched Bifidobacterium longum Strains. npj Biofilms Microbiomes 11, 91 (2025). https://doi.org/10.1038/s41522-025-00714-w
Received:
Accepted:
Published:
DOI: https://doi.org/10.1038/s41522-025-00714-w
