Disordered development of gut microbiome interferes with the establishment of the gut ecosystem during early childhood with atopic dermatitis.

The gut microbiome influences the development of allergic diseases during early childhood. However, there is a lack of comprehensive understanding of microbiome-host crosstalk. Here, we analyzed the influence of gut microbiome dynamics in early childhood on atopic dermatitis (AD) and the potential interactions between host and microbiome that control this homeostasis. We analyzed the gut microbiome in 346 fecal samples (6-36 months; 112 non-AD, 110 mild AD, and 124 moderate to severe AD) from the Longitudinal Cohort for Childhood Origin of Asthma and Allergic Disease birth cohort. The microbiome-host interactions were analyzed in animal and in vitro cell assays. Although the gut microbiome maturated with age in both AD and non-AD groups, its development was disordered in the AD group. Disordered colonization of short-chain fatty acids (SCFA) producers along with age led to abnormal SCFA production and increased IgE levels. A butyrate deficiency and downregulation of GPR109A and PPAR-γ genes were detected in AD-induced mice. Insufficient butyrate decreases the oxygen consumption rate of host cells, which can release oxygen to the gut and perturb the gut microbiome. The disordered gut microbiome development could aggravate balanced microbiome-host interactions, including immune responses during early childhood with AD.

Introduction

The gut microbiome during early life influences immune system development that affects health in later life.[1,2] Gut microbiome perturbations during early infancy increase the risk of developing allergic diseases in children.[3,4] Several studies support the notion that the gut microbiome plays a critical role in manifesting allergic diseases.[5-8] Gut microbiota, and the short-chain fatty acids (SCFA) they produce, can induce Treg cells, which control mucosal Th2 inflammation.[9,10] The Th1/Th2 balance is essential for immune regulation, and an imbalance Th1/Th2 can lead to chronic inflammation and allergic diseases.[11]

Atopic dermatitis (AD) is the most common chronic inflammatory skin disease, and its course is affected by the microbiome.[8,12-15] The gut microbiome in infants with AD has low bacterial diversity, is deficient in Bifidobacterium and Bacteroides,[8,16,17] and has perturbed functional genes related to host immune development.[13] These perturbations result in imbalanced SCFA production, which dysregulates host-microbiome communication through immune and metabolic processes.[18] Previous studies have compared the gut microbiome between infants with and without AD; however, these studies did not consider the changes in the microbiome with aging. Microbes colonize the gut since birth,[2,19] and the microbiome dynamically changes until 36 months of age.[2] The changes in the complex and mutual relationships between the gut microbiome and hosts with AD with aging remain unclear.

A recent study investigated the longitudinal changes in the gut microbiota in children with AD at early (5, 13, 21, and 31 weeks) and school ages (6–11 years) using 16S rRNA amplicon sequences.[20] Although changes in microbiota associated with AD were observed, there was a lack of understanding of the microbiome functions as the study was solely based on 16S rRNA amplicon sequencing.[21] Since the gut microbiome has highly balanced symbiotic interactions with the host,[22] particularly in early life, a comprehensive analysis of the gut microbiome’s role in early childhood is a prerequisite to identify the mechanism underlying the host-microbiome symbiotic ecosystem in AD.

We investigated the relationship between gut microbiome dynamics in early childhood in children with AD and without AD and the potential interactions that control this relationship.

Results

Gut microbiota during early childhood exhibited two age-specific types (based on human study)

We analyzed 346 fecal samples (6–36 months; non-AD = 112, mild AD = 110, and moderate to severe AD = 124) from the Longitudinal Cohort for Childhood Origin of Asthma and Allergic Disease birth cohort (Figure 1a and Figure S1). The sample size had a discriminatory power of 90%, with a significance of 0.05 (Table S1). Of these, 58.0% were boys, and 25.6% were born by cesarean delivery (Table S2). The scoring of atopic dermatitis (SCORAD) index, parental history of allergic diseases, serum eosinophils, total IgE, egg-, and milk-IgE were significantly different between non-AD and AD groups (p < .001). Solid food was introduced at a mean age of 5.6 months in the subjects of this study.

Figure 1.

Summary of the analyzed data and comparison of the gut microbiota between clusters determined by Dirichlet multinomial mixture (DMM) modeling. (a) Overview of the studied samples. (b) Number of samples analyzed for each dataset. We analyzed 346 16S rRNA gene amplicon sequencing (16S), 343 whole metagenome sequencing (WMS), 326 quantitative real-time PCR (RTP), and 343 short-chain fatty acids (SCFA) profiles. (c) Microbiota clusters determined by DMM modeling in WMS and 16S dataset. The lowest Laplace approximation indicated that the number of clusters was two for gut microbiota data (inner graph). Two gut microbiota types (GMT1 and GMT2) were clearly distinguished in the NMDS plot (P < .001). (d) Comparison of age, bacterial amounts, and diversity between GMT1 and GMT2. (e) Heatmap of the relative abundance of the top 20 genera within two clusters. Mod-Sev AD: moderate to severe AD.*P < .05, **P < .01, ***P < .001.

We analyzed 346 16S rRNA amplicon sequences (16S), 343 whole metagenome sequences (WMS), 326 quantitative real-time PCR (RTP) tests, and 343 SCFAs profiles in this study (Figure 1b and Table S3). Moreover, negative controls for the sampling tube, DNA extraction kit, and library preparation were sequenced and subsequently compared to the microbiota in samples (Figure S2). The microbiota detected in the negative controls differed from that in the samples, indicating that our sequence data was not influenced by potential contaminations.

Gut microbiota variation in all subjects (non-AD and AD groups) was analyzed using the Dirichlet multinomial mixture (DMM) model based on WMS and 16S datasets (Figure 1c). The microbiota was clustered within two gut microbiota types, GMT1 and GMT2, based on the lowest Laplace approximation (p < .001). Age was the significant factor distinguishing the two GMTs (Figure 1d; the median age of 10 months in GMT1 and 24 months in GMT2, p < .001), but AD diagnosis was not significantly different between GMTs (p > .05; Table S4). There were no significant differences in the odds ratio (OR) for covariates between GMT1 and GMT2 (p > .05) except age. The relative amounts and diversity of bacteria were higher in GMT2 than GMT1 (p < .01).

Bifidobacterium, Veillonella, and Escherichia were dominant genera in GMT1, whereas Bacteroides, Bifidobacterium, and Faecalibacterium were dominant in GMT2 (q < 0.001; Figure 1e, Table S5). These coincided with the DMM modeling results in each phenotypic group (Figure S3). Although WMS and 16S datasets showed similar results for DMM clustering, the relative abundances of each genus were different. Therefore, we compared the overall concordance of relative abundances between WMS and 16S at phylum and genus levels (Figure S4). Pearson correlation R values of WMS and 16S were over 0.52 in all subjects (p < .001) at the phylum level and decreased (R ≤ 0.31; p < .001) at the genus level. Moreover, detected genera in negative controls for WMS were clearly different from those in samples, and detected numbers were lower (2 genera) than in 16S (19 genera) (Figure S2). Therefore, we primarily used the WMS dataset for further analyses.

Age-dependent assembly of gut microbiota was perturbed in childhood with AD (based on human study)

Although the early childhood gut microbiota was clustered within two types, previous studies reported that the gut microbiota dynamically changes with age during early life.[2,23] Furthermore, solid food intake starting at a mean age of 5.6 months in subjects of this study could influence their gut microbiome (after 6 months). Therefore, we analyzed the gut microbiota using the Bray–Curtis dissimilarity according to age range (6, 7–12, 13–24, 25–36 months) in non-metric multidimensional scaling (NMDS) plots (Figure 2a). The microbiota varied with age within GMT1 and GMT2 for all subjects and each phenotype (p < .01).

Figure 2.

Age-dependent shifts of the gut microbiota in each group. (a) Variation in the gut microbiota according to age in all subjects, non-AD, mild AD, and moderate to severe AD groups were analyzed using NMDS plots based on the species-level Bray–Curtis dissimilarity. The boxplot shows the variation of Bray–Curtis dissimilarity among age groups for each axis. (b) Shifts of indicator species for each age throughout early childhood. Heatmap shows the comparison of the indicator’s relative abundance according to age. Color codes for the age groups and the characteristics of oxygen dependency are listed below the figure. (c) Comparison of the relative abundance of each indicator among groups according to age. Mod-Sev AD: moderate to severe AD. *q < 0.05.

We analyzed the indicator species for each age to study the shifts in the gut microbiota with age in each phenotypic group (Figure 2b, Table S6). Indicator species gradually shifted with age in the non-AD group, while they showed abrupt changes in AD groups. In the non-AD group, facultative anaerobes were indicators until 12 months, and strict anaerobes were indicators after 12 months. However, their relative abundances in the AD groups were different from those in the non-AD group with age. The relative abundance of facultative anaerobes in AD microbiomes was higher than in non-AD microbiomes during early life. In addition, they differed between mild and moderate to severe AD groups (p < .05; Figure S5a).

At 6 months, Bifidobacterium longum was a predominant indicator in all groups, and there were no significantly different species according to AD severity (q > 0.05; Figure 2c). B. bifidum and B. breve were dominant indicators at 7–12 months, and the relative abundance of Escherichia coli was higher in the moderate to severe AD than mild AD group at this age (q < 0.05). The higher proportions of Alistipes finegoldii, unclassified (UC)_Oscillibacter, and UC_Alistipes were detected in the non-AD than in AD groups at 13–24 months (q < 0.05). Faecalibacterium prausnitzii, Bacteroides fragilis, Ruminococcus bromii, and Dialister invisus were indicators at 13–36 months, and the relative abundances of B. fragilis and D. invisus were higher in the non-AD than in AD groups at 25–36 months (q < 0.05). Indicators for each age were evaluated using the multivariate association with linear models (MaAsLin2) after adjusting for covariates (Table S7). Indicators were significantly associated with age after adjustments except for B. bifidum and A. finegoldii.

We analyzed the significantly associated covariates with microbiota at each age using the EnvFit model (Table S8). Breastfeeding and AD severity were significantly associated with gut microbiota at 6–24 months. In particular, exclusive breastfeeding (EBF) was the most significant and common factor (r[2] = 0.306, p < .001 in WMS) for the gut microbiota differences in both non-AD and AD groups at 6 months. The relative abundance of B. longum was higher in EBF infants than in non-EBF infants, whereas B. breve and UC_Coprobacillus were higher in non-EBF infants (p < .05; Figure S5b).

Estimated microbiota age and network analyses revealed disordered gut microbiome development in AD groups (based on human study)

We compared the species profiles by age to analyze the gut microbiome development during early life according to AD severity via random forest modeling. Bacteroides fragilis and F. prausnitzii were the most influential species in the prediction model (Figure 3a). Gut microbiota in the AD group was highly mature compared to that in the non-AD group at 6 months, whereas delayed gut microbiota maturation was observed after 12 months in the microbial-by-age z-score (MAZ) analysis (Figure 3b). Bacteroides fragilis, UC_Subdoligranulum, Anaerostipes hadrus, and Roseburia inulinivorans were significantly correlated to low MAZ in the AD group after 12 months (p < .05; Figure S6a).

Figure 3.

Disordered gut microbiota development in AD groups. (a) Species importance in the prediction model of EMA in 112 non-AD samples. The 16 species with the most discriminating power were selected by the lowest cross-validation error (inner graph). (b) Comparison of EMA through early childhood among groups by MAZ score. (c) Heterogeneity of gut microbiota within each group throughout early life. The gray area represents the 95% confidence intervals (CIs). The black dashed line indicates the time point at which the change in the MAZ score of AD groups occurs. (d) Difference of Shannon diversity index among groups in each age. (e) Network analysis of the top 10 species in each microbiota type for each phenotype determined by random forest modeling. The common group included commonly detected species in GMT1 and GMT2. Positive correlations are marked by blue edges and negative correlations by red edges. Edge thickness denotes FastSpar correlation, ranging from value −0.4 to 0.4. Node sizes were scaled on the eigenvector centrality measure. Only significant correlations with P < .05 are shown. The top 5 network hubs determined by the PageRank algorithm are marked in bold with an Orange square symbol. Mod-Sev AD: moderate to severe AD. *P < .05.

The microbiota heterogeneity was the lowest in the moderate to severe AD group at 6 months (Figure 3c) but increased to the level of other groups after 12 months. Heterogeneity of microbiota in the mild AD group was higher than in other groups after 24 months. The Shannon diversity increased with age in all groups (p < .001), and the diversity was significantly lower in the moderate to severe AD than in non-AD at 6 months (p < .05; Figure 3d). Differences in diversity between groups according to AD severity were observed at 6 months and 25–36 months in the 16S dataset, and 7–12 months and 25–36 months in the functional gene dataset (p < .05; Figure S6b).

Different microbiota development in AD groups could be partly driven by interactions between species in the gut microbiota. Thus, we analyzed the correlations between species in GMT1 and GMT2 of each phenotype group because gut microbiota was clearly distinguished into two age-specific types throughout early childhood (Figure 1). Commonly detected bacteria in both GMTs could be intermediate members during gut microbiome development.

More complex interspecies correlations in the non-AD group than in the AD groups indicated that dynamic drift and diversification occurred through interactions during gut microbiome development in the non-AD group (Figure 3e). Positive correlations were more abundant in GMT1 than in GMT2 and the common groups. Enterococcus faecalis was a hub with positive correlations with bacteria in GMT1. UC_Subdoligranulum and Lachnospiraceae bacterium_5_1_63FAA were keystones during gut microbiome development. Faecalibacterium prausnitzii and B. fragilis were hubs in the gut microbiota of GMT2 with negative correlations with GMT1 members and positive correlations with bacteria in GMT2.

Citrobacter freundii and E. faecalis were keystones in GMT1 of mild and moderate to severe AD, respectively. Although UC_Subdoligranulum and Lachnospiraceae were also keystones during microbiome development in AD groups, there were fewer correlated bacteria than in the non-AD group, and their correlations differed. Faecalibacterium prausnitzii was also a hub in GMT2 of AD groups. Therefore, the same bacteria could be critical to gut microbiome development in all groups. However, their influences on the microbiota were different in AD groups according to severity. This could be related to the disordered gut microbiome development with lower relative abundances of B. fragilis and UC_Subdoligranulum in AD groups (Figure S6a). Network analysis for all subject and combined AD groups revealed similar results as that of the distinguished phenotypes (Figure S7).

Disordered microbiota development in AD groups altered SCFA profiles with age (based on human study)

Differences in AD gut microbiome development can cause alterations in SCFA profiles. Therefore, we compared SCFA profiles among groups by training a machine-learning algorithm. SCFA-by-age z-score (SAZ) was calculated to analyze the dysbiosis of SCFA profiles in AD groups along with age (Figure 4a). Butyrate contributed more significantly to the prediction model than propionate. The SCFA profiles in AD groups were highly matured compared to those in non-AD before 12 months. However, they delayed maturation after 12 months, consistent with gut microbiota development (correlation R value between SAZ and MAZ ≥ 0.77 with p < .001). These results indicated that disordered gut microbiota development was associated with dysregulated SCFA production.

Figure 4.

Dysbiosis of short-chain fatty acids (SCFAs) profiles in AD groups throughout early childhood. (a) Dysbiosis of SCFAs profiles in AD groups according to age. The alteration of SCFAs in AD groups was analyzed using the prediction model of SCFAs changes according to age in 112 non-AD samples. The importance of butyrate and propionate for the prediction model was determined using subjects without AD. The SCFAs-by-age z-score (SAZ) was calculated to identify dysbiosis of SCFA profiles in AD groups along with age. The gray area represents the 95% confidence intervals (CIs). The correlation between SAZ and MAZ was determined by linear regression models and the Pearson correlation test. (b) Species that are significantly associated with the dysbiosis of SCFAs (SAZ score) in AD groups at each age. (c) Comparison of gene families involved in butyrate metabolism among groups. The boxplot shows the significantly different gene families between the non-AD and AD groups. (d) Correlation between the relative abundances of significantly associated species with abnormal MAZ or SAZ in AD groups and IgE levels. Bacteroides fragilis and unclassified (UC)_Subdoligranulum were detected as associated species with P-value < 0.01. Mod-Sev AD: moderate to severe AD. N: non-AD, M: mild AD, and MS: moderate to severe AD. *q < 0.05, **q < 0.01.

The relative abundances of B. breve in AD groups were significantly associated with high SAZ at 6 months (p = .012; Figure 4b). Conversely, the relative abundances of facultative anaerobes (K. pneumoniae, C. freundii, and E. coli), B. fragilis, and F. prausnitzii were related to low SAZ in AD groups after 12 months (p < .05).

Gene families involved in butyrate metabolism were compared according to AD severity (Table S9 and Figure S8), and the significantly different gene families between non-AD and AD groups (q < 0.05 in Dunn’s test and q < 0.25 in MaAsLin2 after adjustments) were analyzed (Figure 4c). L-glutamate usage for butyrate production through glutamate decarboxylase (gad) was lower in AD groups than in non-AD group after 12 months. The conversion between fumarate and succinate by succinate dehydrogenase/fumarate reductase (sdhB/frdB) was lower in moderate to severe AD at 7–12 months. Different genes were used to convert pyruvate to acetyl-CoA in mild AD (porB) compared to other groups (korA) at 25–36 months. Crotonoyl-CoA production from 3-butenoyl-CoA by 4-hydroxybutanoyl-CoA dehydrase (abfD) was reduced in AD groups at 13–24 months. Butyryl-CoA: acetate-CoA transferase (but) levels at 25–36 months were lower in AD groups than non-AD group, and acetaldehyde dehydrogenase (adhE) levels were higher in the mild AD group than in other groups. Phosphate butyryltransferase (ptb) and butyrate kinase (buk) levels were the lowest in the mild AD group at 25–36 months. Ten species (A. findgoldii, A. hadrus, B. fragilis, B. ovatus, Coprococcus eutactus, Eggerthella lenta, Lachnospiraceae bacterium, Megamonas hypermegale, Odoribacter splanchnicus, and Prevotella buccae) were significantly different between non- and AD groups, as corresponding bacteria (q < 0.05; Figure S9). These differences could cause deficient butyrate production and low SAZ in AD groups. Differences according to severity were also found. Gene families involved in propionate metabolism were also different among groups according to severity (Table S10 and Figure S10).

Potential influence of disordered microbiome development on host immune responses in AD groups was revealed by functional features (based on human study)

The associated species with abnormal MAZ and SAZ in AD groups were correlated with clinical features (Figure 4d). The relative abundances of B. fragilis and UC_Subdoligranulum were negatively correlated with total IgE, egg-, and milk-IgE at 25–36 months (p < .01). In particular, B. fragilis was related to both abnormal MAZ and SAZ in AD groups, suggesting that the disordered microbiome development and SCFA production by limited colonization of B. fragilis was related to IgE levels and sensitization to egg or milk in children with AD.

The gut microbiome may play different roles by disordered microbiota development in AD groups. We analyzed the alteration of functional features in the gut microbiome along with age based on KEGG Orthology (KO). The differential KOs by age in the non-AD group were selected by MaAsLin2 after adjustment for covariates. A total of 327 KOs significantly altered with age (q < 0.05). We selected the top 19 most discriminatory features by random forest model with the lowest cross-fold validation error (Figure 5 and Table S11). Differences in functional features were distinguished by two clusters consistent with DMM clustering. Relative abundances of 6 KOs were higher in the GMT1 dominant phase (6–12 months), and those of 13 KOs were higher in the GMT2 dominant phase (13–36 months).

Figure 5.

Normalized abundances of KEGG Orthology (KO) among groups according to age. (a) The top 19 KOs with the most discriminative power were selected. Functional gene features in the gut microbiome of the non-AD group were used to determine significantly changed KOs according to age by MaAsLin2 after adjustment with other covariates (q < 0.05). Among the age-dependent KOs, 19 KOs were selected by the lowest cross-fold validation error. (b) Selected KOs were compared among groups at each age in the heatmap. The colors of the heatmap correspond with the z-score for normalized abundances of functional genes. (c) Comparison of contributing species to significantly different KOs at each age. Statistical significance was calculated using the Kruskal–Wallis test and Dunn’s test. P-value was adjusted using the Benjamini–Hochberg method to calculate the q value. The gradually changing KOs according to AD severity are indicated by significance (P < .05 in the Kruskal test and q < 0.05 in Dunn’s test). Mod-Sev: moderate to severe AD. *q < 0.05, **q < 0.01.

Four KOs gradually changed according to AD severity at each age (p < .05 in Kruskal–Wallis test and q < 0.05 in Dunn’s test). The relative abundance of acarbose and validamycin biosynthesis was higher in non-AD than in AD groups at 6 months. The MAPK signaling pathway showed the highest levels in moderate to severe AD at 7–12 months. Chloroalkane and chloroalkene degradation and transcription were higher in moderate to severe AD at 13–24 months. MAPK signaling pathway showed higher levels in non-AD than in AD groups at 25–36 months. Species contributed to these KOs differed among the AD severity groups.

Interactions between gut microbiome and host cell were perturbed in the AD group (based on animal and in vitro cell studies)

The interactions between the perturbed gut microbiome and host cell were analyzed using animal and in vitro cell assays. The gut microbiota, SCFA production, and host gene expression were compared between non-AD and AD-induced mice (Figure 6a). The concentration of butyrate was lower in fecal samples of the AD group than the non-AD group, whereas acetate was higher in the AD group (p < .01; Figure 6b). The expression level of genes for G-protein coupled receptor 109A (Gpr109a) and peroxisome proliferator-activated receptor-γ (Pparg) were downregulated in the colon of the AD group compared to the non-AD group (p < .01; Figure 6c). Positively correlated gut microbes with Pparg, Gpr109a, and butyrate were more abundant in the non-AD group than in the AD group (p < .05; Figure 6d).

Figure 6.

The interactions between gut microbiome and host cell in animal and in vitro cell assays. (a) Design for the mice experiment. Female BALB/c mice were used, 1% DNCB was applied to the shaved dorsal skin three times for one week. 0.4% DNCB was applied three times a week for four weeks. Mice were sacrificed at day 35, and the middle part of the colon and fecal samples were collected. (b) The concentration of SCFAs was compared between non-AD and AD groups. SCFAs were measured by GC-MS from fecal samples. (c) The expression of target genes in the middle part of the colon was compared between non-AD and AD groups. The expression levels of SCFA receptor (Gpr109a, Gpr43, and Gpr41), transcription factors (Foxp3 and Pparg), anti-inflammatory cytokine (Il10), and pro-inflammatory cytokines (Tnfa, Il1b, and Il6) were analyzed. (d) Significantly correlated genera in fecal samples of mice with Pparg, Gpr109a, acetate, and butyrate. (e) Effects of butyrate on oxygen consumption rate (OCR) of Caco-2 cells. The OCR of each treated cell was investigated using a Seahorse Bioscience XF24 extracellular flux analyzer. Basal respiration (Basal resp.), ATP-linked respiration (ATP prod.), and maximal respiration (Max. resp.) were compared among treated butyrate concentrations. Results are presented as the mean ± SD. *P < .05, **P < .01.

Our results clearly show that butyrate production by the gut microbiome in the AD group was disturbed in both human and animal studies; hence, the influence of butyrate on colon epithelial cells was analyzed by the oxygen consumption rate (OCR) in Caco-2 cells (Figure 6e). Butyrate enhanced OCR and mitochondrial respiration processes. However, this effect was only detected at 1 mM concentration of butyrate, which showed that an optimal concentration of butyrate is required to maintain homeostasis of host cell metabolism.

Discussion

We analyzed the gut microbiome of children with AD during early childhood (6–36 months) according to severity. The ecological drift of the gut microbiome during early life, from facultative anaerobes to strict anaerobes, was irregular in AD groups. The disordered microbiome development in AD was contributed to butyrate deficiency, different microbiome’s functional genes, and increased IgE levels. Perturbed gut microbiota and deficient butyrate interacted with gene expressions and OCR in host cells. The perturbed microbiome-host crosstalk could contribute to AD during early childhood.

The gut microbiota in early childhood was clustered into two age-specific types, GMT1 and GMT2, regardless of AD phenotype. Age-dependent clustering is consistent with previous findings,[2] resulting from synergistic influences, such as immune system development and the transition to solid food during early life.[19,24]

The age-dependent gut microbiome development during early childhood is critical for understanding their influences on AD, since microbial diversity and specific taxa could be linked to AD. Some studies have reported that the gut microbiota of AD patients is less diverse than non-AD patients.[16,25] In contrast, other studies have found no significant differences in microbiota diversity between the groups.[13,26] Bifidobacteria are suggested to be beneficial,[27,28] but are overrepresented in children with allergies at a later age.[7] Faecalibacterium, Lachnospira, Veillonella, and Rothia in the microbiota of children aged 3 months are risk factors for developing atopic wheeze, but not at 12 months.[29] We found that the diversity of gut microbiota in moderate to severe AD decreased compared to the non-AD group at only 6 months through WMS. Bifidobacterium was not associated with AD, as previously reported.[13,20,30,31] Therefore, the disturbed gut ecosystem could be more critical in AD progression during early childhood than single taxa.

Although the gut microbiota developed with age in both non-AD and AD groups, the development in the AD groups was disordered, as determined by the MAZ analysis. We found that the microbiota of children with AD over-matured before 12 months compared to those of non-AD children and had delayed development at 13–36 months. Shifts of indicator species at each age and interspecies interactions indicated the disordered gut microbiome development in AD groups. The order and timing of microbe colonization affect interspecies interactions during gut microbiome development.[32]

Early colonizers can affect the subsequent colonization of gut microbial populations. Interspecies interactions in AD microbiomes were reduced compared to the non-AD microbiome, particularly in GMT1, and these interactions varied according to AD severity. Also, the reduced interactions in GMT1 persisted in GMT2. This provides evidence for the influence of priority effects on the gut microbiome during early development. The low heterogeneity and diversity of moderate to severe AD gut microbiota at 6 months partially supports this assertion. Reduced interspecies interactions at early phases caused reduced interactions at intermediate and later phase (GMT2). This can cause disordered microbiome development and abruptly shift indicator species along ages in the AD groups.

The SCFA dysbiosis with age in AD groups (SAZ analysis) was consistent with MAZ results. Lower SCFA levels and butyrate producers in early life are associated with atopic disease development.[14,33,34] The maturation of the gut microbiome and its metabolome with age are a synchronized process that regulates gut maturation and immune development.[35] Therefore, SCFA dysbiosis resulting from disordered gut microbiome development during early childhood can be related to AD.

The relative abundance of facultative anaerobes throughout early childhood was higher in AD groups than in non-AD group in the present study. Early colonizing facultative anaerobes gradually shift to strict anaerobes by utilizing oxygen to create an anaerobic gut environment.[19,36] The higher abundance of facultative anaerobes in AD groups than in non-AD group could indicate that the redox balance of the gut environment throughout early childhood is disturbed for those with AD. Low SAZ score in AD groups was related to the relative abundances of facultative anaerobes and strict anaerobes including B. fragilis.

Lower B. fragilis proportions were also related to low MAZ scores, reduced interspecies interactions, and increased IgE levels in AD groups compared to non-AD. In particular, IgE and specific IgE to food allergen play a critical role in the pathogenesis of AD and food allergy.[37-39] A reduced abundance of B. fragilis in the AD group, which prevents inflammation through its polysaccharide A by restoring Th1/Th2 balance,[40] was also found in previous studies.[16,27] These results suggest that the redox-dysregulated gut environment of children with AD during early development could cause gut microbiome perturbation and SCFA dysbiosis in later life by limiting the colonization of B. fragilis. This could cause increasing IgE and egg/milk-IgE levels in AD.

The potential influence of disordered gut microbiome development on host immune responses in AD groups was found in functional features analysis. Functional genes in the microbiome related to acarbose and validamycin biosynthesis and MAPK signaling pathway were gradually changed according to AD severity. The anti-inflammatory effect of acarbose has been previously reported in the mouse model of inflammation.[41] MAPK signaling pathway was higher in AD groups at the earlier phase, whereas this pathway was higher in non-AD group at the later phase. This pathway is important for cell growth, differentiation, and inflammation,[42,43] and therefore may play different roles in the gut through different contributing species at early and later phases of early childhood.

The interaction between the perturbed gut microbiome and host cell in the AD group was evaluated in animal studies and in vitro cell assays. Deficient butyrate was caused by perturbed gut microbiota in the AD-induced mice, resulting in the downregulated expression of Gpr109a and Pparg in the colon. GPR109A is a receptor for butyrate,[44] and butyrate activates PPAR-γ signaling, which induces epithelial cell proliferation, oxidative phosphorylation, and immune development.[45-47] In addition, butyrate enhanced the OCR of the host cell, which corresponds to an increase in mitochondrial oxidative phosphorylation and a decrease in anaerobic glycolysis.[48,49] When the host’s metabolism shifts toward anaerobic glycolysis, it decreases oxygen consumption and increases glucose consumption and lactate release,[50,51] driving the proliferation of facultative anaerobes.[52] These influence the differentiation of Treg cells and immune responses.[10,53] This circular crosstalk can accelerate gut ecosystem dysbiosis. Therefore, the disordered gut microbiome development and SCFA dysbiosis in children with AD affect gut homeostasis, including the immune system. Our results advance the understanding of the gut microbiome-host interactions in AD during early childhood (Figure 7).

Figure 7.

Hypothetical model for microbiome-host interactions in the gut ecosystem of AD during early life. In the non-AD group, the gut microbiome develops during early childhood by colonizing microbes at the appropriate time. Early colonizing facultative anaerobes gradually shift to strict anaerobes by utilizing oxygen to create an anaerobic gut environment. Adequate maturation of the gut microbiome regulates SCFAs production and immune development throughout early childhood. Butyrate produced by strict anaerobes activates PPAR-γ, mitochondrial respiration, and increases oxygen consumption through oxidative phosphorylation. Thus, it maintains an anaerobic environment by reducing oxygen emanation from the mucosa. Moreover, the butyrate induces Treg cells, which control Th2 inflammation. In the AD group, this cyclic interaction is imbalanced, exacerbating the homeostasis of the gut ecosystem. An abnormal oxygen environment in the gut maintains facultative anaerobes during childhood; it limits B. fragilis colonization throughout the lifespan. Disordered gut microbiome development is related to decreased butyrate production and abnormal immune responses. PPAR-γ and Treg cells cannot be activated with a decrease in butyrate. The epithelial proliferation and mitochondrial respiration decrease, and the emanating oxygen and lactate from mucosa increase by anaerobic glycolysis of the host cell. Facultative anaerobes use lactate produced by the host cell. Abnormal Treg and epithelial cell proliferation and the increased numbers of facultative anaerobes induce abnormal immune responses.

Since it is challenging to evaluate human gut microbiome development using animal models,[54,55] we could not apply microbiome development states in mice. This can be improved with advancements in omics technologies, in vitro, and humanized in vivo models in the future. Although differences of microbiome’s functional gene among phenotypic groups were found by whole metagenome in this study, further studies will be necessary to clarify these findings by other techniques including metatranscriptomic and metaproteomic analyses. Furthermore, we could not identify the causative covariates of gut microbiome development in this study. This could be due to the limited covariates obtained or complex interactions in the gut ecosystem. Nevertheless, our study is significant because it advances the hypothesis regarding the influence of the gut microbiome on the pathogenesis of AD through host-microbiome interactions during early life. This can help develop novel prediction and treatment strategies for AD in early life.

In conclusion, regulated gut microbiome development is essential to maintain the gut ecosystem during early childhood. Disordered microbiome development in AD is characterized by persisting facultative anaerobes and limited B. fragilis colonization with age, which reduces SCFA production and induces abnormal immune responses by increasing IgE. This relationship is dynamic and harmonious crosstalk of a symbiotic human-microbiome system. Thus, the early stages of microbiome development may be an appropriate target for modulating the microbiota and establishing a healthy microbiome. Understanding these symbiotic relationships is critical for developing preventive and therapeutic strategies to treat AD in early childhood.

Material and methods

Study subjects and sample collection

The statistical power of the sample number per group was estimated using the micropower R-package based on permutational multivariate analysis of variance (PERMANOVA).[56] The PERMANOVA powers were calculated based on the weighted UniFrac distance of infant gut microbiota data in our previous study.[13] The effect sizes (ω[2]) were calculated using the simulated matrixes of 80% and 90% powers for varying sample numbers per group (Table S1). One hundred bootstrap iterations were performed using a significance level of 0.05 to estimate the power and effect size. The effect size was smaller than 0.0001 in ≥ 100 samples per group for a discriminatory power of 90%. Therefore, we aimed to collect more than 100 samples per group in this study.

Subjects in this study were recruited from the Longitudinal Cohort for Childhood Origin of Asthma and Allergic Disease (COCOA) birth cohort (over 3,000 infants enrolled) and Childhood Asthma Atopy Asan Medical Center.[57] The collection of data and biological samples in the COCOA study is conducted every year for all enrolled children from 6 months to 20 years old, regardless of allergic disease development. Pediatric allergists diagnosed AD according to Hanifin and Rajka’s criteria.[58] AD severity was simultaneously assessed at fecal collection using the SCORAD index (mild < 25 and moderate to severe ≥ 25).[59] Total IgE and egg- and milk-specific serum IgE (IU/mL) levels were measured using fluorescent enzyme immunoassay (AutoCAP system, Phadia AB, Uppsala, Sweden) after 12 months of age for all subjects. Subjects without any visible sign of skin eczema indicative of AD from 6 to 36 months and without food sensitization (< 0.35 of specific IgE to egg and milk) were recruited into the non-AD group. Subjects who received antibiotics during 6 months preceding collection and had health complications, including gestational age < 37 weeks, smoking exposure during pregnancy, any major congenital anomalies, and birth asphyxia requiring oxygen supplementation, were excluded from the study. Collected clinical data and biological samples in this study were summarized in the sampling strategy (Figure S1).

A total of 346 samples (112 non-AD, 110 mild AD, and 124 moderate to severe AD) from 6 to 36 months old were collected from the COCOA cohort. Fecal samples were collected from five university hospitals in Seoul, and each center transported samples in iceboxes to the laboratory within 4 hours after sample collection. Samples were immediately stored at −80°C before being processed for DNA extraction.

This study was approved by the Institutional Review Boards (IRBs) of Asan Medical Center (IRB no. 2008–0616 and 2015–1031), Samsung Medical Center (IRB no. 2009–02-021), Severance Hospital (IRB no. 4–2008-0588), CHA Medical Center (IRB no. 2010–010), and Seoul National University Hospital (IRB no. H-1401-086-550). Written informed consent was obtained from the parents of each infant.

16S rRNA gene amplicon sequencing

Metagenomic DNA was extracted from fecal samples using the RNeasy PowerMicrobiome Kit (Cat #26000-50, Qiagen, Valencia, CA, USA) following the manufacturer’s instructions. The extracted DNA was purified using the DNeasy PowerClean Pro Cleanup Kit (Cat #12997-50, Qiagen) and quantified using a BioPhotometer D30 with a μCuvette G1.0 (Eppendorf, Hamburg, Germany). The V1–V3 region of the 16S rRNA gene was amplified using a C1000 thermal cycler (Bio-Rad, Hercules, CA, USA) per the MiSeq system protocol for preparing a 16S metagenomics sequencing library (Illumina, Inc., San Diego, CA, USA), as described previously.[60,61] Equimolar concentrations of each sample were pooled and sequenced using the Illumina MiSeq system (300-bp paired ends) according to the manufacturer’s instructions. Previous studies reported the influence of potential contamination of reagents on sequence-based microbiome analyses.[62,63] Therefore, we sequenced and analyzed negative controls for quality control of sequencing. Negative controls were included at every step to check contamination, and three negative controls were sequenced with samples. Sequenced negative controls included unused stool box (sampling blank), DNA-free water added to the RNeasy PowerMicrobiome Kit (negative extraction controls), and the library preparation instead of DNA (library negative controls).

Amplicon sequences were analyzed using the QIIME2 pipeline.[64] Raw sequences were quality filtered, and denoised using DADA2, and the taxonomic position of representative sequences were assigned with the EzTaxon-e database.[65] Diversity indices were calculated after rarefied without replacement. A total of 28,804,166 reads (median 68,045 reads per sample) in human fecal samples were obtained from sequence analyses.

Whole metagenome shotgun sequencing

Extracted DNA from 346 samples was fragmented using a NEBNext dsDNA Fragmentase (Cat #0348 L, New England Biolabs, Ipswich, MA, USA). Then, metagenomic libraries were prepared using the ACCEL-NGS 2S PLUS DNA Library Kits (Cat #21096, Swift Biosciences, Ann Arbor, MI, USA) according to the manufacturer’s instructions.

The size of the libraries was confirmed using a Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA). The concentration of the library was measured using a PicoGreen dsDNA Assay kit (Cat #P11496, Invitrogen, Carlsbad, CA, USA). Equimolar concentrations of each library were calculated by qPCR using a TaKaRa PCR Thermal Cycler Dice Real Time System III (TaKaRa Bio, Inc., Shiga, Japan) with the GenNext NGS Library Quantification Kit (Cat #NLQ-101, Toyobo, Osaka, Japan).

Libraries were pooled and sequenced using the Illumina HiSeq 2500 system (250-bp paired ends). Three samples were excluded due to failed sequencing library preparation. Thus, a total of 343 samples were analyzed. For quality control of shotgun sequencing, negative controls were processed following the same procedures to check contamination, and two kinds of negative control were sequenced with samples. The two negative controls were DNA-free water added to the RNeasy PowerMicrobiome Kit (negative extraction controls) and the library preparation instead of DNA (library negative controls).

WMS were analyzed as described previously.[13,66] Briefly, adapter removal and quality filtering were performed using Trimmomatic with default options (leading and trailing: 3, sliding window size: 4, quality: 15, min length: 150).[67] Paired-end sequences were merged using PEAR v.0.9.11.[68] Contaminated human genes were removed using the BBMap (http://sourceforge.net/project/bbmap) with a reference human genome. Taxonomic features were obtained by using MetaPhlAn2 v.2.7,[69] and functional features were obtained by using HUMAnN2 v.2.8.0.[70] Finally, the resultant UniRef90 IDs were converted to the KOs. A total of 2,635,653,581 reads (median 6,837,236 reads per sample) were obtained from WMS analysis.

Quantification of total bacterial amounts

The relative amounts of bacteria in fecal samples were estimated by quantitative real-time PCR based on the 16S rRNA gene. The 16S rRNA gene was amplified with the primer 340 F (5′-TCC TAC GGG AGG CAG CAG-3′) and 518 R (5′-ATT ACC GCG GCT GCT GG-3′) using a Thermal Cycle Dice Real-Time System III (TaKaRa Bio, Inc.).

Each sample was measured in triplicates in a 25-μL reaction containing 12.5 μL of TB Green Premix Ex Taq (Tli RNaseH Plus) (Cat #RR820B, TaKaRa Bio, Inc.), 2 μM of each primer, and 1 μL of DNA template (a 10-fold dilution series of sample DNA) with following amplification conditions: 95°C for 30s, followed by 40 cycles of denaturation at 95°C for 5 s and annealing at 60°C for 30s. We quantified the bacterial amount by comparing threshold cycles (Ct) values to a standard curve generated from parallel reactions of serial dilutions (1 × 101–1 × 107) of the 16S rRNA gene from the Escherichia coli K12 w3110 strain. Regression coefficients (r[2]) for all standard curves were ≥ 0.98.

SCFAs analysis

Fecal SCFAs were extracted with modification of previously described methods.[71,72] Acetic acid (Cat #338826, Sigma-Aldrich, St. Louis, MO, USA), propionic acid (Cat #81910, Sigma), iso-butyric acid (Cat #58360, Sigma), butyric acid (Cat #W222100, Sigma), iso-valeric acid (Cat #129542, Sigma), and valeric acid (Cat # 240370, Sigma) were used as standard compounds. 4-methlyvaleric acid (Cat #277827, Sigma) was used as an internal standard (IS) for final adjustment on the quantity of each compound in samples. Briefly, 20 mg of fecal sample was suspended in extraction solution (450 μL of methylene chloride, 5 μL of IS (10 ppm), and 25 μL of 0.6 M HCl). The suspension was derivatized for 10 min after vortexing and centrifuging (14,000 × g) for 10 min at 4°C. Finally, the methylene chloride layer was transferred into the test vial and briefly vortexed before injecting into the gas chromatography-mass spectrometry (GC-MS).

The GC-MS analysis was performed on an Agilent 7890A gas chromatograph coupled to an Agilent 7000 triple quad mass spectrometric, detector (MSD, Cat #G7000-90038, Agilent Technologies). Derivatives were separated using a VF-WAXms capillary column (30 m × 0.25 mm, 0.25 μm film thickness, Agilent J & W Scientific, Folsom, CA, USA) with a carrier gas (helium) at a flow rate of 2.5 mL/min. One microliter of the sample was injected into the system. The oven program was set at 80°C for 2 min, raised to 160°C at a rate of 20°C/min, to 180°C at a rate of 3°C/min, to a final temperature of 230°C at a rate of 10°C/min and maintained for 5 min. The transfer line temperature was set to 200°C. The MS source was held at 230°C and the quadrupole at 150°C. The energy of electron ionization was set to 70 eV. The mass spectral data were collected in a selected ion monitoring mode.

Quantitative analysis software (MSD, Agilent Technologies) was used to process the GC-MS data for peak picking, standard curve construction, and SCFAs quantification. The concentration of each SCFA in samples was calculated using the calibration curve constructed from the GC-MS data of corresponding SCFA standards.

Predictive SCFAs metabolite profiling

We used the Model-based Genomically Informed High-dimensional Predictor of Microbial Community Metabolic Profiles (MelonnPan v.0.99.0)[73] to predict SCFAs profiles based on measured SCFAs in fecal samples and functional gene features of gut microbiomes obtained from HUMAnN2 (UniRef90 gene families).

After normalizing the raw measures into relative abundances, we filtered features (gene families and SCFAs) by both relative abundance (> 0.01%) and prevalence (> 10% of the samples). We measured SCFAs and functional gene features from 41 samples to construct a training set. SCFAs were predicted based on the functional gene features of 302 samples using elastic net regularization and cross-validation. The predictability of each metabolite was evaluated using the Spearman correlation coefficient (r) between the measured and predicted metabolites concentration across samples.

The same labels were repeatedly shuffled in both metabolite and gene family tables to test the significance of well-predicted metabolites (r > 0.3),[73] the MelonnPan model was applied to the randomized data by shuffling to link genes to metabolites, and the predicted metabolites in the original data and randomized data were compared.

Butyrate and propionate were successfully predicted after evaluating 10-fold cross-validation of predicted metabolisms in the present study. The accuracy of predicted metabolites was evaluated by the representative training sample index (RTSI) scores based on principal component analysis (PCA).

Animal studies

Female C57BL/6 mice were purchased from ORIENT Bio Korea. The mice (8-week-old) were divided into AD (n = 8) and non-AD (n = 10) groups. AD was induced for five weeks. Briefly, 200 μL of 1% 1-Chloro-2,4-dinitrobenzene (DNCB) in an acetone:olive oil mixture (3:1 vol/vol) was applied to the shaved dorsal skin, and then the skin was covered with a transparent film dressing. 1% DNCB was applied three times a week for one week. Then, 0.4% DNCB was applied three times a week for four weeks. At 35 days, mice were sacrificed, and the middle part of the colon was collected for gene expression analysis, and fecal samples were collected for microbiota and SCFA analyses. All procedures were performed in accordance with the guideline of the Institutional Animal Care and Usage Committee at Soonchunhyang University.

Quantitative reverse transcriptase-polymerase chain reaction analysis

Total RNA was extracted using TRIzol reagent (Cat #15596018, Invitrogen) from the collected middle part of the colon. RNA was converted to cDNA using reverse transcription reagents (Cat #FSQ-201, Toyobo). PCR was performed using SYBR Green Real time PCR Master Mix Kit (Toyobo) with specific primers for target genes (Table S12). The reaction was performed with QuantStudio5 Real-Time PCR System (Applied Biosystems, Foster City, CA, USA). The expression levels of target genes were calculated by comparing the relative expression levels after normalization to β-actin.

Microbiota analysis in mice fecal samples

Metagenomic DNA was extracted from fecal samples of mice using the RNeasy PowerMicrobiome Kit (Qiagen). The 16S rRNA gene amplicon sequencing and analyses were performed as described above. A total of 442,006 reads (median 20,543 reads per sample) in mouse fecal samples were obtained from sequence analyses.

SCFAs analysis

SCFAs in mice fecal samples were extracted and analyzed using GC-MS as described previously.

In vitro cell experiment

Caco-2 cells (ATCC® CRL-1790TM) were cultured in Eagle’s minimum essential medium (EMEM) supplemented with 20% fetal bovine serum (FBS) and 1% penicillin-streptomycin by incubation at 37°C in a humidified atmosphere containing 5% CO2. Cultured Caco-2 cells at passages 20–27 were used for experiments.

Cell viability

Caco-2 cells were seeded (0.5 ×104 cells/well) into 96-well plates and incubated at 37°C for 24 h to allow the cells to attach. Caco-2 cells were treated with butyrate (concentrations of 0.5, 1, and 5 mM). Afterward, 20 μL of MTT solution (final concentration of 1 mg/mL) was added to each well, and cells were incubated for 2 h. The cell culture medium was subsequently removed, and each well was treated with 200 μL DMSO to dissolve formazan crystals. Dissolved formazan absorbance was measured at 570 nm using a microplate reader (Sunrise-Basic Tecan, Tecan Austria GmbH, Grödig, Austria). The cell viability was expressed as the percentage of MTT reduction calculated relative to the absorbance of control cells.

Mitochondrial Function

The effect of butyrate on mitochondrial function was measured by seeding Caco-2 cells in Seahorse XFp mini-plate (Cat #103025-100, Agilent Technologies) at 0.7 ×104 cells/well. Cells were cultured as described above. The cells were washed with extracellular flux (XF) DMEM medium (Agilent Technologies) supplemented with 5.55 mM glucose, 2 mM glutamine, and 1 mM sodium pyruvate followed by incubation in this medium for 60 min at 37°C in a non-CO2 incubator. Plates were transferred to a Seahorse XFp analyzer (Agilent Technologies) and subjected to an equilibration period. After measuring the basal oxygen consumption rate (OCR) for four cycles, oligomycin (1.5 µM, Agilent Technologies) was added to determine the proportion of respiration used to generate ATP. Then, carbonyl cyanide-4-(trifluoro methoxy) phenyl hydrazone (FCCP, 0.5 µM, Agilent Technologies) was added to determine the maximal respiration by mitochondria. After that, rotenone (0.5 µM, Agilent Technologies) and antimycin A (0.5 µM, Agilent Technologies) were added to measure the non-mitochondrial respiratory rate.

Statistical analysis

Statistical analysis was performed with the R software v.4.0.2. All statistical tests for microbiome data were two-sided.

Clinical data

The Kruskal–Wallis test was used to compare continuous variables among groups, and the Chi-square test was used to compare categorical variables between groups. The ORs and corresponding 95% confidence intervals (CI) were calculated using unconditional logistic regression models.

Microbiota analysis

The difference in Shannon diversity index between non-AD and AD groups was analyzed using the Mann–Whitney–Wilcox test. Differences in beta-diversity were visualized using NMDS plots and tested for inference by PERMANOVA (Adonis from the package vegan with 999 permutations) based on the Bray–Curtis distance.

Differences in relative abundances in microbiota were analyzed using the Kruskal–Wallis tests for three-group comparisons, and Dunn’s test was used to identify pairwise differences between groups. Dunn’s test was performed using the ‘dunn.test’ package in R, and an approximately normal distribution was used to calculate p-values. P-values were adjusted using Benjamini–Hochberg false discovery rate (FDR) multiple testing, correction. Results with q (adjusted p-value) < 0.05 were considered statistically significant. WMS data was normalized using the cumulative sum scaling (CSS) method. Differences in normalized abundance by CSS among groups were analyzed using multivariate association with linear models (MaAsLin2).[74]

The inter-individual gut microbiota heterogeneity within the same group according to age was determined by the locally weighted scatterplot smoothing (LOESS) regression based on the Bray–Curtis dissimilarity. The heterogeneity was evaluated using the Kruskal–Wallis test and Mann–Whitney U test. The significance of quantitative real-time PCR results among groups was calculated using the Mann–Whitney U test and Kruskal–Wallis test.

The correlations between significantly associated genera with Pparg, Gpr109a, acetate, and butyrate in mice experiments were analyzed using the Spearman correlation in R software.

DMM clustering

The variance of the gut microbiota in all samples was analyzed using DMM modeling[75] with the R package DirichletMultinomial. DMM clustering was conducted at the genus level to compare results between 16S and WMS. The lowest Laplace approximation score determined the number of clusters.

Comparison of taxonomic profiles between 16S and WMS

The overall concordance of microbiota compositions between 16S and WMS analyses was compared at the phylum and genus levels. Relative abundances of taxa were visualized using bar plots. Linear regression models and a Pearson correlation test were used to investigate the coefficient inference associations between 16S and WMS datasets. The correlation and regression of two datasets were visualized using scatterplots.

Indicator species analysis

The indicator species at each age were determined using the ‘multipatt’ function of the indicspecies package in R software. We selected species with > 10% prevalence and > 0.01% mean relative abundance in at least one group to identify indicators. This procedure is a statistical method determining whether a particular species is significantly more abundant in predefined groups than when the same species are randomly assigned to the groups. The non-AD microbiota determined the indicator species at each age. The significance of the association of the indicator with age was tested using 10,000 permutations. Results with p < .05 were considered indicators at each age. Resultant indicators were evaluated by multivariate analysis of MaAsLin2 after adjustment for covariates.

EnvFit analysis

Covariates (sex, delivery mode, feeding type, allergic family history, AD diagnosis, total IgE, egg-IgE, and milk-IgE) could be related to differences in the gut microbiome and allergic disease. The effect size and significance of each covariate in the variation of microbiota were determined using the ‘envfit’ function in R package vegan (v.2.5–7), which compared the difference in centroids of each group relative to the total variation. The significance was determined using 999 permutations. A PERMANOVA was used to determine the covariates with the strongest effects on taxonomic and functional features. Results with p < .05 were considered statistically significant.

Random forest model

Random forest regression was performed to model the gut microbiota age based on the relative abundance of bacterial species obtained from WMS of 112 non-AD samples (6–36 months) using the R package Ranger. The chronological age of the microbiota in children without AD was used as a training set to predict the estimated microbiota age (EMA). The resulting prediction model, defined by alterations in bacterial species, was subsequently applied to all subjects with Ranger’s ‘predict’ function. These estimates were used as a proxy for gut microbiota maturation.[76] Next, sensitivity analyses were performed by restricting the models to samples not included in model building to confirm the independence of the training sets.

To estimate the microbiota age, the N top discriminatory species were identified by the ‘rfcv’ function in the randomForest package. We calculated the MAZ to compare the maturation of the gut microbiota among groups as described previously.[20,76] A lower MAZ is indicative of a delay in microbiota development.

We used a similar process to estimate the SAZ for analyzing the dysbiosis of SCFAs profiles according to gut microbiota development variations between non-AD and AD groups, and predicted the SCFA age score using the training set. Linear regression models and a Pearson correlation test were used to analyze the associations between the SAZ and MAZ. In addition, Spearman correlation was used to identify the species associated with SAZ in each phenotypic group.

The correlations between significantly associated species with abnormal MAZ or SAZ in AD groups and clinical features (SCORAD index, Total IgE, egg-IgE, milk-IgE, and eosinophil) were analyzed using the Spearman correlation in R software.

Network analysis

Interspecies correlations were estimated using FastSpar[77] based on Pearson correlation with 1,000 bootstraps. Pseudo p-values were calculated as the proportion of simulated bootstrapped datasets with a correlation at least as extreme as the one computed for the original data set. The corresponding correlation network was visualized using the R package qgraph. The top 10 species in each microbiota (GMT) group were selected for each phenotype by random forest, and they were used for network analysis. Significant correlations (with p < .05) were shown in the network.

Node sizes were scaled on the eigenvector centrality measure, which was determined using the R package igraph. Edge thickness denoted a FastSpar correlation ranging from values −0.4 to 0.4. Hubs were identified using the PageRank algorithm, a link analysis with the underlying assumption that hubs are more connected to other nodes than non-hub nodes.[78] The top five species with the highest PageRank in each phenotypic group were selected as hubs in the networks. This process was performed in R using the igraph package.

Functional features analysis

Functional features of the gut microbiome were analyzed using the KO category. The significantly different features, according to age, were selected by MaAsLin2 after adjustment for covariates (q < 0.05). The N top discriminatory KOs were selected by ‘rfcv’ function in the randomForest package.

Gene families involved in butyrate and propionate metabolism were compared between phenotypic groups at each age. MaAsLin2 evaluated the difference in gene families between phenotypic groups at each age after adjusting for covariates. We evaluated the significance of different KOs and gene families between phenotypic groups at each age using the Kruskal–Wallis test and Dunn’s test. P-values were adjusted using the Benjamini–Hochberg method. Results with q < 0.05 were considered statistically significant.

Click here for additional data file.