Changes in Gut Microbiota of Patients with Atopic Dermatitis During Balneotherapy

Introduction

Atopic dermatitis (AD) is a common inflammatory skin disorder affecting 15–30% of children and 2–10% of adults.¹ AD results from a complex interaction between genetics and the environment both decreasing skin barrier function. Abnormalities of the skin barrier, the immune barrier, and more recently, the microbial barrier or microbiota have been described.²

Thermal Spring Water (TSW) has been used successfully since the eighteenth century in treating skin diseases including atopic dermatitis (AD) and psoriasis.^3,4 La Roche-Posay thermal spring water (LRP-TSW) in France is a selenium-rich water (0.053 mg/l) that contains bicarbonate, calcium, silicate, magnesium, and strontium.⁵ The regimen consists of an 18-day treatment (3 weeks) including a daily high-pressure filiform shower (15 bars for 3 minutes) performed by a dermatologist. Clinical trials have demonstrated positive effects of balneotherapy with LRP-TSW on AD.^5,6 Recently, the role of LRP-TSW to improve the diversity of the skin microbiota and to reduce dryness and pruritus in AD and psoriasis has been highlighted.⁶ The concentration of minerals and non-pathogenic microbes in LRP-TSW may explain its therapeutic benefit.³ Clinical studies have shown that topical LRP-TSW treatment results in increase in Gram-negative bacteria with reduction of Gram-positive bacteria, and improvement in skin microbial diversity.⁷ At the same time, skin condition and general dryness have been shown to improve.⁸

Gut microbiota was particularly investigated during the last decade for its potential role in immune-mediated dermatological diseases such as AD through the gut-skin axis.^9–12 This cross-talk between gut and skin could be mediated by immunity, by metabolites (such as short-chain-fatty acids) or by neuroendocrine molecules produced by the gut microbiota.¹⁰ In attempting to characterize the gut microbiota of patients with AD, observational studies, considering especially infants or children, used various techniques (mainly 16S rRNA sequencing, but also qPCR or cultivation) and reported sometimes-contradictory results.¹³ Gut microbiota diversity was either decreased in patients with AD,^14–16 increased,¹⁷ or without any difference.¹³ Most studies focused mainly on Bifidobacteria, Clostridia and Lactobacilli.¹³ Depending on the studies, Bifidobacteria were either increased^13,18 or decreased^10,13,19 in AD patients. The same contradictory results were observed with Clostridia,^16,20,21 while Lactobacilli were mostly found decreased.^13,18 Other studies reported that AD gut microbiome showed higher abundance of Clostridioides difficile, Escherichia coli and Staphylococcus aureus.¹⁰

In this study, we used shotgun metagenomics to evaluate gut microbiota changes during BT, and how it is connected with the improvement of the pathology observed after BT (SCORAD).

Materials and Methods

Patient Recruitment and Sampling

This open label observational study did not need an institutional review board approval but complied with the ethical guidelines of the 1975 Declaration of Helsinki and was conducted between April and September 2016 according to ICH guidelines for Good Clinical Practice. This research used anonymous data without direct patient contact, and institutional review board approval was not necessary prior to the study initiation. Ninety-six individuals, 51 females and 45 males, between 6 and 70 years old, suffering from AD were followed before and at the end of an 18 days balneotherapy treatment at the thermal care center of La Roche-Posay (France) (Table 1). Patients were not prospectively assigned balneotherapy but rather as part of routine treatment, as deemed by the dermatologists. All subjects provided informed consent to participate in the study, when patients were minors, consents were signed by their parents/guardians instead. The electronic file compiling the data recorded in the study was submitted to the French Commission in charge of Personal Data Protection and was approved for data treatment. The study protocol and corresponding appendices were also transmitted to a French Research Ethics Committee (CPP, Poitiers, France) for information purposes. The CPP confirmed the non-necessity to get an approval for this study.

Table 1 Phenotypic Description of the Cohort and Its Evolution During the Balneotherapy

During the 2 visits (D0 and D18), the same investigating dermatologist assessed the severity of AD with the SCORAD (SCORing Atopic Dermatitis) based on SCORAD A (disease extension), SCORAD B (disease intensity) and SCORAD C (subjective symptoms such as pruritus or sleep loss).²²

At the same visits, fecal samples were collected in a stabilizing solution (RNAlater^®) according to the International Human Microbiome Standards (IHMS), SOP 005 V2.²³

DNA Extraction and Shotgun Metagenomics Sequencing

Fecal DNA was extracted following the SOP 07 V2 H from IHMS.²³ The DNA preparation was subjected to quality control using Qubit Fluorometric (ThermoFisher Scientific, Waltham, US) and qualified using DNA size profiling on a Fragment Analyzer instrument (Agilent Technologies, Santa Clara, US).

3 µg of high molecular weight DNA (>10 kbp) was used to build sequencing libraries. Shearing of DNA into fragments of approximately 150 bp was performed using an ultrasonicator (Covaris, Woburn, US) and DNA fragment library construction was performed using the Ion Plus Fragment Library and Ion Xpress Barcode Adapters Kits (ThermoFisher Scientific, Waltham, US). Purified and amplified DNA fragment libraries were sequenced using the Ion Proton Sequencer (ThermoFisher Scientific, Waltham, US) leading to 22.7 million reads (± 1.3 million) on average per sample.

Sequencing Data Processing

Read cleaning, filtering and mapping were performed with the METEOR software suite²⁴ that relies on Bowtie2 for read mapping.²⁵ First, quality control was performed with AlienTrimmer:²⁶ sequencing adapters were removed and low quality reads were trimmed or discarded. Then, reads mapped to the human genome (GRCh38.p13) with nucleotide identity ≥90% were also discarded.

Remaining reads were mapped to the updated Integrated Gene Catalogue of the human gut microbiota (IGC2),²⁷ comprising 10.4 million of genes. Gene counts were computed with a two-step procedure as previously described.²⁸ The resulting count table was further processed using the R package MetaOMineR v1.31.²⁹ Downsizing at 15 million mapped reads was performed to take into account differences in sequencing depth and in mapping rate across samples. Then, the downsized matrix was normalized according to gene length and transformed into a frequency matrix (FPKM normalization).

MetaGenomic Species (MSP) Abundance Profiles

The IGC2 gene catalogue has been previously clustered into 1990 MetaGenomic Species (MSP, clusters of >100 co-abundant genes belonging to the same microbial species, in 1601 samples) using MSPminer.³⁰ Taxonomic annotation of MSP was performed using an in-house pipeline. First, all genes are aligned on public databases (NCBI nt and WGS, versions of March 2020³¹) using blastn.³² An MSP was assigned to the lowest taxonomic rank (from species to superkingdom) that brought consensus in at least 50% of its genes. To avoid misleading annotations due to error in databases, for each gene the 20 first hits were being considered. Relative abundance of an MSP was computed as the mean abundance of its 100 “marker” genes (that is, the genes that correlate the most altogether). If less than 10% of “marker” genes were seen in a sample, the abundance of the MSP was set to 0. Relative abundances at higher taxonomic ranks were computed as the sum of the MSP that belong to a given taxa. MSP richness was assessed as the number of MSP detected in a sample (that is, whose abundance is strictly positive).

Functional Profiling

Three databases were used to predict gene functions: Kyoto Encyclopedia of Genes and Genomes (KEGG);³³ eggNOG database³⁴ and TIGRFAM.³⁵ Genes from the IGC2 catalogue were mapped with diamond³⁶ onto KEGG Orthologous (KO) from KEGG database (version 89). Each gene was assigned to the best-ranked KO among hits with e-value < 10^e-05 and bit score ≥60. The same procedure was used with eggnog (version 3.0). TIGRFAM families (version 15.0) were searched in the IGC2 catalogue using HMMER 3.2.1.³⁷

Then we assessed presence of KEGG modules, Gut-Metabolic Modules³⁸ and Gut-Brain Modules³⁹ in an MSP. A functional module consists in a set of KOs (or NOGs, or TIGRFAMs). Since MSPs are pangenomes, their genes are divided into “core” genes which are present in all samples harboring the MSP or “accessory” genes which might be absent from a sample even if the MSP is detected. Thus, we first considered a functional module to be present in an MSP if at least 90% of its components were present in the “core” genes of the MSP. Then we re-affined this assumption sample by sample, by adding to the “core” genes the accessory genes detected in a given sample. Finally, we estimated the abundance of a module in a sample by summing abundances of all MSPs found to carry this module in the sample.

Statistical Analysis

All statistical analyses and graphs were performed with R software (version 3.6.0).⁴⁰ For each individual, we computed the SCORAD reduction between D18 an D0 as:

With being the SCORAD at D18 and the SCORAD at D0.

We defined as richness species’ MSPs that were positively correlated to richness, and ‘anti-richness’ species MSPs that were negatively correlated to richness (at baseline, p ≤ 0.05, Spearman’s correlation). Differential analysis in MSPs or functional modules abundances were performed using Wilcoxon signed-rank tests. Correlations between variables (either metagenomics variables or clinical variables) were performed using Spearman’s correlations. All p-values were corrected for multiple testing with the Benjamini Hochberg method. Unless stated otherwise, a corrected p-value is considered as significant when under the threshold of 0.1.

To account for age effect in the relations between disease duration and metagenomic features, we computed Spearman’s partial correlations using the R package ppcor v1.1.⁴¹

Principal Components Analysis (PCA) were performed on MSPs relative abundance after they have been log10-transformed with the package FactoMineR v1.42.⁴² Clinical variables were used as supplementary variables.

Effect size was computed using the Cliff’s Delta statistic with the package effsize v0.7.4.⁴³ This metrics compares abundance of a given value between two groups and is comprised between −1 and 1 (0: no difference between both groups regarding the considered value; from 0 to 1: the value is larger in the first group as compared to the second group; from −1 to 0: the value is larger in the second group as compared to the first group). Magnitude of the effect size d is defined in⁴⁴ as negligible if |d| < 0.147, small if |d| < 0.33, medium if |d| < 0.474, and large otherwise.

Results

Abundances of Species Associated with High or Low MSP Richness Correlate with SCORAD and SCORAD Reduction

At D0, disease severity estimated by SCORAD was not directly associated with the number of metagenomic species detected in samples (MSP richness). Indeed, correlation between MSP richness and SCORAD was not significant (Spearman’s ρ = −0.022, P = 0.82, Figure S1). However, species associated with high richness (“richness” species) tended to correlate negatively to SCORAD and positively to SCORAD reduction, whereas species associated with low richness (“anti-richness” species) tended to correlate positively to SCORAD and negatively to SCORAD reduction (Figure 1A-C). Among the “anti-richness” species, we found Escherichia coli, Streptococcus sanguinis, two Eggertellaceae (Eggerthella lenta and Rubneribacter badeniensis), two Lachnospiraceae (Sellimonas intestinalis and Anaerotignum lactatifermentans), and two Erysipelotrichaceae (Longicatena caecimuris and Faecalitalea sp.) that were all positively correlated to SCORAD (P ≤ 0.05, Spearman correlation), and Anaerotruncus colihominis that was negatively correlated to SCORAD reduction (P ≤ 0.05, Spearman correlation). In contrast, a group of “richness” species, composed of six Lachnospiraceae, three Ruminococcaceae (among which Faecalibacterium prausnitzii, subspecies msp_0457), one Eggerthellaceae (Senegalimassilia anaerobia), one Christensenellaceae and two unclassified Clostridiales, correlated negatively with SCORAD (P ≤ 0.05, Spearman correlation). In addition, another group of “richness” species correlated positively with SCORAD reduction (among which two Lachnospiraceae, two Eggerthellaceae, two Christensenellaceae) (P ≤ 0.05, Spearman correlation).

Figure 1 (A) Spearman correlations between abundance of Metagenomics Species (MSP) at D0 and MSP richness, SCORAD (D0) and SCORAD reduction. Only MSPs with at least one significant correlation with SCORAD or SCORAD reduction, and significantly correlated with gut richness, are shown. Black dots represent correlations with FDR ≤ 0.1; empty circles represent correlations with p ≤ 0.05. (B and C) Relation between MSP correlations with richness and SCORAD at baseline or SCORAD reduction; (D) Relation between species correlation with SCORAD at baseline and reduction. Spearman correlations are indicated.

MSPs that tended to be negatively associated with disease severity also consistently tended to be positively associated with SCORAD reduction, and vice versa (Figure 1D). Two species had significant associations: an unclassified Lachnospiraceae and an unclassified Christensenellaceae, both “richness” species associated with lower SCORAD at baseline and with SCORAD reduction.

Moderate but Positive Changes in Gut Microbiota During Balneotherapy

At baseline, mean SCORAD of the 96 individuals was 65 ±13, with 84 individuals (87.5%) suffering from severe AD (SCORAD ≥ 50), and 12 individuals (12.5%) suffering from moderate AD (25 ≤ SCORAD < 50) (Table 1). The 18-day balneotherapy drastically reduced the mean SCORAD (by 33.8%), going down from 65 ±13 to 43 ± 14 (P < 2.2e-16, Wilcoxon signed-rank test, Table 1, Figure 2A). Among the 84 individuals with severe AD at D0, 55 (65%) had moderate AD at D18, and 6 (7%) had slight AD. 33% of the 12 individuals with moderate AD at D0 had slight AD at the end of BT and the remaining 67% had still moderate AD. At the end of the treatment, there were only 23 individuals (24.0%) with severe AD while 63 (65.6%) had moderate AD and 10 (10.4%) had mild AD (SCORAD < 25). However, we found no difference between MSP richness at baseline and end of treatment (p = 0.51, Figure 2B).

Figure 2 (A) SCORAD evolution during the balneotherapy. (B) MSP richness evolution during the balneotherapy. (C) Cliff’s Delta of contrasted MSPs (above) and functional modules (below) between D0 and D18, along with their correlation to SCORAD reduction (Cliff’s Delta positive: enriched at D0; Cliff’s Delta negative: enriched at D18). (D) Correlations between MSPs abundance (at D0) and SCORAD amelioration (y-axis) in function of the Cliff’s Delta between D0 and D18 of the MSPs (x-axis). (E) Correlations between modules abundance (at D0) and SCORAD amelioration (y-axis) in function of the Cliff’s Delta between D0 and D18 of the modules (x-axis). Empty circles represent correlations with p ≤ 0.05.

At the species level, differences between D0 and D18 were slight. Out of 554 MSP present in at least 10% of individuals, 23 (4%) were differentially abundant between D0 and D18 (P ≤ 0.05, Figures 2C and S2). Among them, 14 were reduced at D18 and 9 were increased in relative abundance. Only Mediterraneibacter faecis, Lactobacillus delbrueckii and Eisenbergiella tayi had a non-negligible effect size (absolute Cliff’s Delta > 0.1), and the 3 of them were depleted at D18. Of interest, Faecalibacterium prausnitzii (phylogroup represented by the L2-6 strain) was enriched at D18, even though the effect size was negligible (Cliff’s Delta = −0.044).

Regarding functional potential of species, 24 out of 270 functional modules (9%) were differentially abundant between D0 and D18 (P < 0.05) (Figure 2C). Noteworthy, 14% of Gut-Brain Modules (GBM), which focus on metabolites involved in the gut-brain axis, were found to be differentially abundant, while only 7% and 8% of GMM and KEGG modules (that also contain modules related to gut-brain axis but are not specific to them) were contrasted, respectively. GABA degradation was increased at D18, while acetate synthesis, quinolinic acid degradation, tryptophan degradation, quinolinic acid synthesis, menaquinone synthesis and nitric oxide synthesis were decreased at D18. Many modules related to tryptophan metabolism were more abundant at D0: tryptophan degradation (MGB049 and MF0009), quinolinic acid synthesis, quinolinic acid degradation, NAD biosynthesis from aspartate, and all of them had Cliff’s Delta >0.1.

The changes in abundance of an MSP between D0 and D18 (measured by Cliff’s Delta) was associated with its relation to SCORAD reduction (Spearman’s ρ = −0.12, P = 0.0042, Figure 2D): MSPs more abundant at baseline in well-improving patients tended to increase during the treatment, while MSPs more abundant in poorly-improving patients tended to decrease. The same global effect was observed at the functional level (Spearman’s ρ = −0.23, P = 0.00013, Figure 2E).

Microbiota Composition is Associated with History of Atopy

MSP richness at D0 and D18 was positively correlated with disease duration (Spearman’s ρ = 0.24, P = 0.017 and ρ = 0.28, P = 0.005, respectively, Figure 3A and B), which means that long-term patients had a richer gut microbiota than short-term patients (Table 1). After accounting for age with Spearman’s partial correlation, the association between disease duration and MSP richness was still significant at D18 (Spearman’s ρ = 0.23, P = 0.02) and close to be at D0 (Spearman’s ρ = 0.18, P = 0.07), suggesting that there was a positive association between disease duration and MSP richness independently of age. This observation was confirmed by Principal Component Analysis (PCA), since disease duration and age projected on orthogonal components (component 1 and component 2, respectively, see Figure S3).

Figure 3 (A and B) MSP richness as a function of disease duration at D0 and D18. (C) Spearman correlations between abundance of Metagenomics Species (MSP) or functional modules, disease duration and SCORAD reduction. Only MSPs and functional modules correlated with disease duration at both D0 and D18 (p ≤ 0.05) are shown. Correlation between disease duration and metagenomics features are computed using a partial Spearman correlation to take age into account. Black dots represent correlations with FDR ≤ 0.1; empty circles represent correlations with p ≤ 0.05. (D) Correlations of MSPs with SCORAD amelioration (y-axis) in function of correlations of MSPs with disease duration (x-axis). (E) Correlations of functional modules with SCORAD amelioration (y-axis) in function of correlations of functional modules with disease duration (x-axis).

The relation between richness and disease duration was of special interest since we observed that species which positively correlated with severity were negatively correlated with richness, and vice versa (c.f. previous section). Inconsistently, we found that disease duration was positively associated with severity at baseline, especially regarding SCORAD B (sub-score of SCORAD related to intensity of symptoms) (Spearman’s ρ = 0.26, P = 0.01, Figure S4). This relationship disappeared at D18 (Spearman’s ρ = 0.09, P = 0.37), suggesting that the long-term patients might be better responders to the treatment. However, disease duration was not significantly correlated to SCORAD reduction (Spearman’s ρ= 0.13, P = 0.20).

At the species level, 33 MSPs were correlated with disease duration at both D0 and D18 after adjusting for age (P ≤ 0.05, Spearman’s partial correlation). Among them, 22 were increased in abundance in long-term patients and 11 were increased in short-term patients (Figure 3C). The species increased in long-term patients tended to be positively correlated to SCORAD reduction while those of short-term patients were negatively correlated to SCORAD reduction (Spearman’s ρ = 0.38, P < 2.2e-16, Figure 3D). The former species included Lachnospira pectinoschiza and Faecalibacterium prausnitzii (phylogroup represented by the strain CNCM I-4573), while the latter species included Clostridium symbiosum, Anaerotignum lactatifermentans, Enterocloster bolteae, Eggerthella lenta and Christensenella minuta. The same pattern was observed at the functional level (Figure 3E). We found 17 modules correlated with disease duration when taking age into account (P ≤ 0.05, Spearman’s partial correlation), among which only one was increased in long-term patients. This unbalanced ratio may be an annotation bias, since 9 out of 11 short-term MSPs (81%) were annotated at the species level vs only 14 out of 22 long-term MSPs (64%). Nevertheless, among the 16 modules increased in short-term patients we found numerous functions related to gut-brain axis, notably glutamine degradation (ammonia pathway), glutamate degradation (4-aminobutanoate pathway), GABA synthesis, menaquinone production, serotonin synthesis and acetylcholine synthesis. All of them were negatively associated with SCORAD reduction (Figure 3C).

We divided the cohort into two subgroups: patients with history of AD ≥ 23 years (n = 48, “long” patients) and patients with history of AD < 23 years (n = 48, “recent” patients). Both subgroups had 42 species whose abundance changed between D0 and D18 (P ≤ 0.1, Wilcoxon signed-rank test). However, only two species were common to long and recent patients, and among them, only one had a consistent evolution between the groups: Desulfovibrio piger, a sulfate-reducing bacterium which prevalence lies below 20% in both groups. Moreover, median Cliff’s delta of the species from the long patients group was −0.02, meaning that most species increased during balneotherapy, whereas median Cliff’s delta of the species from the recent group was +0.02, meaning that most species decreased during balneotherapy. This difference in median was significant (P = 0.03, Wilcoxon test, Figure S5).

Discussion

In this study, we have shown that species associated with high gut microbiota richness were negatively associated with SCORAD and positively associated with SCORAD reduction. Since dysbiotic gut microbiota has been linked to a low richness and inflammatory markers,²⁹ these results suggest that the inflammatory potential of the gut microbiota, AD severity, and symptoms reduction are connected altogether. In particular, the “richness” species Faecalibacterium prausnitzii (subspecies msp_0457), a butyrate-producing bacterium with documented anti-inflammatory properties^45,46 was negatively associated with SCORAD and positively associated with SCORAD reduction. On the opposite, Escherichia coli, “anti-richness” species with LPS associated pro-inflammatory properties and here associated with disease severity, was found to be more abundant in infants with AD in several studies.¹⁰ However, no direct correlation of gut microbiota richness and baseline SCORAD was detected, possibly because of the insufficient statistical power of the study.

Changes in the gut microbiota during the 18-days balneotherapy were slight, with only 4% of the gut species exhibiting a significant but moderate or negligible evolution (according to the effect size). However, abundance of these species at baseline was related to SCORAD evolution, in a consistent way: increasing species were positively associated with SCORAD reduction, while decreasing species were negatively associated. In other words, increasing species were more abundant at baseline in individuals showing the highest SCORAD reduction, while decreasing species were more abundant at baseline in individuals showing the lowest SCORAD reduction. These results suggest that gut microbiota might co-evolve with disease severity. Of interest, a subspecies of Faecalibacterium prausnitzii (msp_0388) increased during the balneotherapy and was positively associated with SCORAD reduction. However, this particular phylogroup of F. prausnitzii was previously found to be more abundant in individuals with AD in a study using 16S gene sequencing.⁴⁷ This difference in results may arise from technical differences (16S vs shotgun). In our study, relative abundance of this species was negatively related to severity (SCORAD and SCORAD B), sleep problem, and itching (data not shown).

One of the decreasing species with the largest effect size (Cliff’s Delta = 0.12) was Lactobacillus delbrueckii, possibly due to a decrease in dairy product consumption during BT. Its baseline abundance was also negatively correlated with SCORAD reduction (Spearman’s ρ = −0.13, P = 0.2). This result is inconsistent with different studies^48–50 that found that two subspecies of L. delbrueckii (bulgaricus and lactis) reduced inflammation in AD, and thus are thought to alleviate AD.

At the functional level, GABA degradation was positively associated with SCORAD reduction and increased during the treatment. On the other hand, modules related to tryptophan metabolism decreased during balneotherapy and were negatively associated with SCORAD reduction. Interestingly, GABA and tryptophan were both found to be related to AD.^51–55

Finally, even though the changes in gut microbiota were slight, there was a significant relation between these changes and the disease evolution, suggesting that gut microbiota and AD tended to co-evolve during the treatment. The short duration of the balneotherapy (18 days) and the absence of medication known to strongly affect the gut microbiota such as antibiotics might explain why these changes were subtle and why effect size was low.

Patients with a long history of atopy were more severe at baseline than patients with a recent disease onset (especially regarding SCORAD B), but this relationship disappeared after balneotherapy. On the other hand, the former exhibited a gut microbiota with beneficial markers such as higher MSP richness and species positively related to SCORAD reduction, while the latter exhibited lower MSP richness and potentially detrimental species. Faecalibacterium prausnitzii (phylogroup represented by the strain CNCM I-4573) was among the species enriched in patients with a long history of atopy, and is known for its anti-inflammatory properties. Enterocloster bolteae, Eggerthella lenta and Anaerotignum lactatifermentans, increased in patients with a recent disease onset, were of particular interest. Indeed, both E. bolteae and E. lenta were found increased in asthmatic individuals.⁵⁶ As for A. lactatifermentans, as well as E. lenta, they were part of the “anti-richness” species which abundance positively correlated with disease severity at baseline. Moreover, gut microbiota displayed a different evolution pattern considering these two groups: most species impacted by the balneotherapy in long-term patients tended to increase at D18, while they mostly decreased in short-term patients. Considering all these observations together – difference in microbiota at baseline, difference in microbiota evolution, disappearance of the baseline acute severity in the long-term patients – we might hypothesize that AD course may induce change in habits,⁵⁷ which in turn affect gut microbiota as well as balneotherapy benefits.

Not only gut-skin axis, but also gut-brain axis might play a key role in the relation between gut microbiota and AD, patients with AD being at increased risk of psychological wellbeing issues.⁵⁸ In our study, we found that 29% of the contrasted functional modules between D0 and D18 were specific to the gut-brain axis modules, whereas they represent only 18% of all considered modules. Among the different modules associated to gut-brain axis, modules related to GABA are of particular interest. Indeed, GABA degradation (MGB019) was significantly associated with SCORAD reduction (positively) and significantly increased between D0 and D18. Moreover, GABA biosynthesis from putrescine (M00135) was significantly associated with disease severity (positively) (with SCORAD A at D0 and at D18, data not shown). Finally, GABA synthesis (MGB022) was more abundant in short-term patients, and negatively associated with SCORAD reduction (though non-significantly). These different results suggest that the ability of the microbiota to degrade GABA could be beneficial for the host, while conversely its ability to produce GABA might have noxious effects on the host in the context of atopy. Since GABA or GABA receptor agonists administration have been shown to alleviate induced-AD in mouse models,^54,55 these findings could seem inconsistent and more GABA-specific studies should be carried out to explore thoroughly the relationship between GABA dosage in the different host’s compartments, GABA intake and gut-microbial GABA production.

Recent studies have shown that skin exposure to external elements, such as ultraviolet⁵⁹ or topical medication⁶⁰ could impact gut microbiota composition. In our study, we also reported changes in gut microbiota occurring during balneotherapy. However, since balneotherapy takes place at La Roche Posay center, it comes along with changes in environmental factors that could confound the observed modifications. For example, diet is known to play a key role in shaping gut microbiota. Absence of information about changes in diet during balneotherapy – including LRP-TSW drinking – is a major limitation in our study. Other uncontrolled factor such as physical exercise or stress reduction might also account for changes in gut microbiota.

Another major limitation is the absence of control group, which was not feasible given the way of recruitment (AD patients coming to the centre to follow the balneotherapy). Hence, the observed changes in patients gut microbiota cannot be bounded to balneotherapy alone but may rather be due to confounding factors such as the mixture of treatment and all environmental factors.

Conclusion

In a balneotherapy treatment for patients with atopic dermatitis, disease parameters such as severity, severity reduction or disease duration, are intricately linked to gut microbiota composition. However, we cannot exclude that changes in the gut microbiota might be due to an effect of lifestyle and diet modifications that occurred during the balneotherapy. Nevertheless, the gut-brain-skin axis, via neurotransmitter such as GABA, should be further studied in diseases such as AD.