Dysbiosis of skin microbiota with increased fungal diversity is associated with severity of disease in atopic dermatitis.

BACKGROUND: Atopic dermatitis (AD) is a multifactorial inflammatory skin disease and an altered skin microbiota with an increase of Staphylococcus aureus has been reported. However, the role of fungi remains poorly investigated.
OBJECTIVES: We aimed to improve the understanding of the fungal skin microbiota, the mycobiota, in AD in relation to the bacterial colonization.
METHODS: Skin swabs of 16 AD patients and 16 healthy controls (HC) from four different skin sites, that is antecubital crease, dorsal neck, glabella and vertex from multiple time points were analysed by DNA sequencing of the internal transcribed spacer region 1 (ITS1) and 16S rRNA gene for fungi and bacteria, respectively.
RESULTS: Malassezia spp. were the predominant fungi in all subjects but with a decreased dominance in severe AD patients in favour of non-Malassezia fungi, for example Candida spp. For bacteria, a decrease of Cutibacterium spp. in AD patients in favour of Staphylococcus spp., particularly S. aureus, was observed. Further, both bacterial and fungal community compositions of severe AD patients significantly differed from mild-to-moderate AD patients and HC with the latter two having overall similar microbiota showing some distinctions in bacterial communities.
CONCLUSIONS: We conclude that severe AD is associated with a pronounced dysbiosis of the microbiota with increased fungal diversity. Potentially infectious agents, for example Staphylococcus and Candida, were increased in severe AD.

Introduction

Atopic dermatitis (AD) is a chronic relapsing inflammatory skin disease with a high impact on quality of life. , While in children the disease mostly affects the creases, in adults head and neck AD and chronic eczema of the hands occur besides the classical pattern. , , The pathogenesis of AD is complex, multifactorial and intertwining. On the one hand, genetic mutations such as the loss‐of‐function mutation in the FLG gene coding for filaggrin, a key regulator of epidermal homeostasis, play a role. On the other hand, impaired skin barrier functions result in an elevated pH and trans‐dermal water loss. Further, the skin immune system is altered in patients with AD shifting towards the type 2 immune response with elevated levels of T helper (Th) type 2 cytokines such as interleukin‐4 (IL‐4), IL‐5 or IL‐13. , ,

Additionally, dysbiosis of the microbes residing on skin is associated with the disease. While healthy skin microbiota composition depend on the skin location, for example with Cutibacterium spp. dominating lipid‐rich areas, skin communities of AD patients are generally reorganized with increased abundance of Staphylococcus aureus and decreased overall microbial diversity. , , , , However, the role of fungal communities in AD remains poorly investigated. The importance of fungal skin microbiota, the mycobiota, is highlighted by the fact that AD patients are frequently sensitized to Malassezia spp., the most abundant fungus on human skin. In contrast to bacteria, fungal diversity appears to increase in AD patients. , To date, no next‐generation sequencing (NGS) studies investigating the mycobiota in patients with different severities of AD have been conducted.

In this study, we aim to improve the understanding of the mycobiota in relation to the bacterial communities by analysing samples of typical lesion sites at multiple time points in AD patients and healthy controls (HC).

Materials and methods

Participant recruitment and sampling

The study was approved by the Ethics Committee of Zurich, Switzerland (BASEC‐no. 2016‐00301). AD patients and HC were adults (>18 years) and all provided written informed consent. AD patients (n = 16) underwent physical examination by a dermatologist at the Department of Dermatology, University Hospital Zurich, Switzerland. The eczema locations of each patient are depicted in Table S1. The severity of disease was determined with the SCORAD system. Averages of SCORAD scores of the different time points were calculated and defined as mild‐to‐moderate (<50) and severe (>50). HC (n = 16) were included if they (i) showed no history of AD or any other skin disease, (ii) had no other chronic medical condition or treatment, (iii) performed negative in skin prick tests and (iv) neither used antibiotics nor antifungals within 6 months prior to sampling. Detailed demographic information was obtained by questionnaires and is summarized in Table S2.

To observe the course of disease, AD patients were sampled at three time points (t1 = week 0, t2 = week 2, t3 = week 4). Samples of HC were taken at t1 and t3. Mean time spans were 14.4 ± 1.9 days between t1 and t2 and 30.1 ± 4.9 days between t1 and t3. At each visit, recent medication was recorded and the skin condition was examined. Subjects were allowed to proceed with their regular skin care habits to guarantee natural skin conditions. Fourteen out of sixteen AD patients applied topical steroids at the sampling day (Table S2) and none of the participants used antibiotics or antifungals within 6 months prior to sampling. The four skin sites antecubital crease, dorsal neck, glabella and vertex were chosen to cover some of the most common eczema sites in adult AD patients. Sampling was performed with flocked swabs (Floqswabs or eSwabs, COPAN, Brescia, Italy) soaked in sterile NaCl (0.9%, Braun, Sempach, Switzerland) by rubbing 4–8 cm2 at each skin site multiple times. Swabs were either frozen in Liquid Amies medium provided with eSwabs and stored at −80°C or directly processed by DNA extraction. Swabs with no skin contact and blank water samples were added as negative controls to detect contaminants.

DNA extraction

Swabs for fungal DNA extraction were added to extraction buffer (1 M of Tris–HCl (pH 8), 50 mmol/L of EDTA (Thermo Fisher Scientific, Rheinach, Switzerland), 0.5% Tween 20 (Bio‐Rad Laboratories, Cressier, Switzerland)) with proteinase K (Roche Diagnostics, Rotkreuz, Switzerland) and predigested at 56°C overnight. Swabs for bacterial DNA extraction were incubated in lysozyme digestion buffer (20 mg/mL lysozyme (Sigma‐Aldrich, Buchs, Switzerland), 25 mmol/L of Tris–HCl (pH 8), 2.5 mmol/L of EDTA (Thermo Fisher Scientific), 1% Triton X‐100 (Sigma‐Aldrich)) at 37°C for 60 min. After mechanical cell wall disruption with 0.5 mm beads (Qiagen, Hilden, Germany) in a Tissuelyser (Qiagen), samples were either processed with Masterpure Yeast DNA Purification kit (Epicentre, LuBioScience GmbH, Zurich, Switzerland) for fungal DNA or Invitrogen Purelink Genomic DNA Mini kit (Thermo Fisher Scientific) for bacterial DNA according to the manufacturers' protocols.

PCR amplification, library preparation and sequencing

For fungal analysis, the ITS1 region was amplified with the primers 18S‐F (5′‐GTAAAAGTCGTAACAAGGTTTC‐3′) and 5.8S‐1R (5′‐GTTCAAAGAYTCGATGATTCAC‐3′) with following cycling conditions: preincubation at 95°C for 300 s, 33 cycles of 98°C for 20 s, 61°C for 20 s, 72°C for 40 s and final extension at 72°C for 60 s. For bacterial analysis, the V1‐3 region within the 16S ribosomal RNA (rRNA) gene was amplified with primers 27F (5′‐AGAGTTTGATCMTGGCTCAG‐3′) and 534R (5′‐ATTACCGCGGCTGCTGG‐3′) to retrieve an accurate composition of the bacterial skin microbiome. PCR conditions were preincubation at 95°C for 180 s, 25 cycles of 95°C for 30 s, 62°C for 30 s, 72°C for 30 s and final extension at 72°C for 300 s. For these first PCRs, 3 μL of extracted DNA was used with the Kapa Hifi Hotstart polymerase (Roche Diagnostics). PCR success was controlled on a 1.5% agarose gel and if necessary, PCR was repeated with 1 and 5 μL of extracted DNA. Libraries were generated according to the instructions of Illumina's 16S rRNA metagenomics sequencing library preparations with some adaptations. Briefly, Nextera XT Indices (Illumina, San Diego, CA, USA) were added by a second PCR with eight cycles. The amplicons of the second PCR were purified with Agencourt Ampure XP beads (Beckman Coulter, Krefeld, Germany), quantified with a Qubit 2.0 Fluorometer (Thermo Fisher Scientific), normalized to 4 nmol/L, and pooled. Further, libraries were prepared with the Miseq Reagent kit v3 (Illumina) and 2 × 300 bp paired‐end reads were generated with the Miseq sequencing machine (Illumina).

Taxonomic classification

Internal transcribed spacer raw sequencing data were processed using PIPITS (v2.3) with default parameters and the RDP Classifier (v2.1211) against UNITE database (04.02.2020) to retrieve taxonomic annotation for each sequence. The taxonomy table was manually adjusted by using BLAST for interesting sequences with >1000 contigs and >97% identity.

16S raw sequencing data were processed using dada2. Briefly, reads were quality filtered (expected errors set to 2 for forward and 3 for backward read), error rates were estimated and reads were merged. Chimeric sequences were removed using the consensus method of the function removeBimeraDenovo. Taxonomic assignment for each ASV was performed using IdTaxa (DECIPHER package) , with RDP v18 as reference database. ASVs not belonging to the kingdom Bacteria or with unassigned phylum were removed and the phylogenetic tree was reconstructed using a general‐time reversible model with GAMMA correction (phangorn package). ASVs with >5000 contigs and >99% identity were manually aligned with BLAST to distinguish between Staphylococcus species.

Data analysis and statistics

The R software was used for data analysis and statistics. For data management of ITS and 16S rRNA data and for plotting the following packages were used: phyloseq, ape, tidyverse including ggplot2 with some add‐ons, gridExtra, magrittr and data.table. Decontamination was performed with the package decontam with the prevalence method and thresholds of 0.1 for fungi and 0.5 for bacteria. For normalization of compositional data, we used corncob (Count Regression for Correlated Observations with the Beta‐binomial) as it has been shown recently that with the frequently used rarefying method valuable information is omitted. , Additionally, corncob was used for calculations based on relative abundances with adjusted P‐values with a False Discovery Rate (FDR) of 0.05. To compare relative abundances of non‐Malassezia fungi in HC and AD, we performed Kruskal–Wallis tests with posthoc Wilcoxon signed‐rank tests. Correlations were performed with Spearman analysis. For α‐diversity, Shannon diversity indices were estimated with DivNet and for statistical downstream analyses Breakaway was used. , To investigate ß‐diversities the PhILR (Phylogenetic Isometric Log‐Ratio Transform) method was used. The data were transformed into an orthogonal unconstrained space related on the phylogenetic tree, which allows graphical display of distances in a principal coordination analysis (PCoA) plot. Permutational multivariate analysis of variance (PERMANOVA) tests of ß‐diversities were performed using the adonis function as implemented in vegan. ,

Results

Study population and samples

We collected 320 skin swabs from 16 AD patients (nine with mild‐to‐moderate and seven with severe disease) and 16 HC from the antecubital crease, dorsal neck, glabella and vertex at two (HC) and three time points (AD). The ITS and 16S rRNA regions were sequenced to investigate the fungal and bacterial communities, respectively. Characteristics of the study population with demographic information, disease background and medication are illustrated in Table S2. All AD patients suffered from typical skin lesions at least at one of the four sampling sites at one or more visits.

Severe atopic dermatitis patients are more frequently colonized with non‐ fungi

First, we investigated relative abundances of fungal communities. The single most abundant fungal class was Malasseziomycetes at all skin sites (Figs 1a, S1a). In severe AD it was significantly decreased compared to HC and mild‐to‐moderate AD (Fig. 2a; Wald test, FDR adjusted P < 0.05). This reduction was in favour of higher relative abundances of non‐Malassezia fungi in severe AD (Kruskal–Wallis and posthoc Wilcoxon tests, P < 0.05) at sebaceous skin sites. Next to Malasseziomycetes, Saccharomycetes including Candida and Debaryomyces was the most abundant class (Fig. 2b).

Figure 1

Relative abundances of fungal and bacterial classes in HC and AD patients. Relative abundances of fungal (a) and bacterial (b) classes with frequencies of >1% merged for healthy controls (HC), mild‐to‐moderate atopic dermatitis (AD) and severe AD at the four different skin sites antecubital crease, dorsal neck, glabella and vertex.

Figure 2

Relative abundances of individual fungal and bacterial taxa in HC and AD patients. Relative abundances of the fungal classes Malasseziomycetes (a) and Saccharomycetes (b), the fungal species M. restricta (c), M. globosa (d), M. sympodialis (e) and M. furfur (f), the bacterial families Staphylococcaceae (g) and Propionibacteriaceae (h), and the bacterial species S. aureus (i) and S. epidermidis (j) grouped for healthy subjects (HC), mild‐to‐moderate atopic dermatitis (AD) and severe AD and the skin sites antecubital crease, dorsal neck, glabella and vertex. Note: The y‐axis was adjusted for less abundant taxa resulting in omission of single data points. Significant differences with false discovery rates of 0.05 are labelled with stars: *P < 0.05; **P < 0.01; ***P < 0.001.

We then analysed the 10 most abundant fungal genera using effect size plots (Fig. S2). In severe AD, the relative abundance of Malassezia was overall decreased, while it was increased at several skin sites for Candida, Debaryomyces, Aureobasidium and Penicillium when comparing to both mild‐to‐moderate AD and HC (Wald test, FDR adjusted P < 0.05).

We further investigated Malassezia at the species level. M. restricta and M. globosa showed the highest relative abundances in all three groups, yet M. restricta was significantly less abundant in severe AD patients (Figs 2c,d, S1b, Wald test, FDR adjusted P < 0.05). The third most abundant species in HC was M. sympodialis that occurred significantly less abundant in severe AD. Vice versa, M. furfur was increased in severe AD being the third abundant species in this group (Figs 2e,f, S1b). Other Malassezia species such as M. slooffiae and M. dermatis dominated the mycobiota of single individuals but were not detected in other subjects (Fig. S1b).

Bacterial dysbiosis is driven by staphylococci and cutibacteria

Next, we studied the bacterial skin microbiota (Figs 1b, S3a). At the antecubital crease, Bacilli including staphylococci was the most abundant class. At the three sebaceous skin sites, Actinobacteria including cutibacteria predominated in HC, while this class was significantly less abundant in AD patients, especially in severe AD (Wald test, FDR adjusted P < 0.05). In contrast, relative abundances of Bacilli were significantly higher in AD patients at sebaceous skin sites, again more pronounced in severe AD (Wald test, FDR adjusted P < 0.05). At the family level, we made likewise observations; a higher relative abundance of Staphylococcaceae and a lower abundance of Propionibacteriaceae in patients with severe AD (Fig. 2g,h; Wald test, FDR adjusted P < 0.05). The three genera Cutibacterium, Staphylococcus and Corynebacterium varied most between AD patients and HC (Fig. S4; Wald test, FDR adjusted P < 0.05). While relative abundance of Cutibacterium was decreased in AD patients, Staphylococcus and Corynebacterium were more abundant.

Finally, we analysed the distribution of Staphylococcus at the species level (Figs 2i,j and S3b). Relative abundances of the two prominent species S. aureus and S. epidermidis were increased in patients with severe AD (Wald test, FDR adjusted P < 0.05). Because the increase of S. aureus in AD patients has been frequently discussed, we further analysed its correlation to the SCORAD (SCORing Atopic Dermatitis). Relative abundance of S. aureus is positively correlated to the SCORAD (Fig. S5).

The fungal diversity increases in severe atopic dermatitis

Shannon diversity indices as a measure of α‐diversity were calculated at the class and genus level in order to analyse both the number of taxa and the inequality between abundances. The fungal diversity was significantly elevated in patients with severe AD compared to HC and mild‐to‐moderate AD at all skin sites except the antecubital crease, where it also was higher in severe AD though only significant in comparison to mild‐to‐moderate AD (Fig. 3a,b). In contrast, the overall bacterial diversity was comparable in all three groups only with a significant increase at the dorsal neck and glabella in mild‐to‐moderate AD (Fig. 3c,d).

Figure 3

α‐Diversities of fungi and bacteria in HC and AD patients. α‐Diversities calculated as Shannon diversity indices (SDI) of fungi (a, b) and bacteria (c, d) at the class (a, c) and genus (b, d) level and grouped for healthy subjects (HC), mild‐to‐moderate atopic dermatitis (AD) and severe AD and the skin sites antecubital crease, dorsal neck, glabella and vertex. Significant differences of pairwise comparison are labelled with stars; *P < 0.05; **P < 0.01, ***P < 0.001.

Microbiota in healthy controls and atopic dermatitis are significantly different

We then investigated the β‐diversities to elucidate the differences of fungal and bacterial community compositions between the three groups (Fig. 4). Overall, the microbiota compositions in HC and severe AD clustered significantly different with the microbiota of mild‐to‐moderate AD being between them. Pairwise statistics performed with PERMANOVA are displayed in Table S3.

Figure 4

ß‐diversities of fungi and bacteria in HC and AD patients displayed in principal coordinates analysis (PCoA) plots. PCoA plots of weighted Euclidean distances in Phylogenetic Isometric Log‐Ratio Transform (PhILR) space for fungal (a) and bacterial communities (b) and grouped for healthy subjects (HC), mild‐to‐moderate atopic dermatitis (AD), and severe AD. Ellipses show the 95% confidence level for a multivariate t‐distribution. Statistics were calculated with pairwise permutational multivariate analysis of variance (PERMANOVA) and are listed in Table S3.

We were further interested whether the compositions were different in lesional vs. unaffected (i.e. non‐lesional) skin in AD patients (Figs S6/S7 and Table S4). However, neither fungal nor bacterial compositions varied between different skin conditions.

Finally, we sampled HC and AD patients at different time points (t1–t3) to elucidate putative short‐term microbial dynamics of the disease. Both fungal and bacterial community compositions did not differ at any skin site over time (Figs S8/S9, Table S5).

Conclusions

In this study, we investigated the skin myco‐ and microbiota with next generation sequencing according to different severities of AD and in HC. We generally confirmed previous studies showing that Malassezia species dominate healthy skin with highest relative abundances of M. restricta and M. globosa. , , , , , Although in severe AD patients Malassezia spp. were also predominant, non‐Malassezia fungi, for example Candida or Debaryomyces, were significantly more abundant than in HC or in mild‐to‐moderate AD. Interestingly, in culture‐based studies, higher frequencies of C. albicans colonization were also found in the oral cavity and faeces of AD patients as compared to healthy controls. , , A significant decrease of M. restricta in severe AD has also been found previously. , However, we are the first to demonstrate a decrease of M. sympodialis and increase of M. furfur in AD. Generally, it is difficult to compare microbiota studies for several reasons: (i) most previous studies did not compare relative abundances of Malassezia species, (ii) sequencing techniques may differ, that is marker gene analysis vs. whole metagenome sequencing, (iii) different targets are used, that is ITS1 vs. ITS2 region, (iv) sampling techniques differ, that is tape stripping vs. swabs, (v) geographical and ethnic differences can lead to altered skin microbiota. , , , , , Thus, to get the whole picture of microbes residing on human skin conducting studies on different continents and standardization of techniques is crucial.

Increased relative abundance of Staphylococcus spp., particularly S. aureus and Staphylococcus epidermidis, as well as decreased Cutibacterium spp. in patients with AD and particularly severe AD were in line with previously published work. , , , , , This dysbiosis was here and elsewhere confirmed by the positive correlation of S. aureus abundance and SCORAD. , , , , The increased abundance of S. aureus might be favourable for the proliferation of Candida spp. as it has been shown in biofilm experiments that these two microorganisms exhibit a synergistic activity. ,

Similar to previous studies, we found an increased fungal diversity in AD patients. , Interestingly, it was only increased in patients with severe AD and not in mild‐to‐moderate cases. We were unable, however, to observe the previously described decrease of bacterial diversity in severe AD. , The bacterial diversity was comparable in all three study groups and was even slightly elevated in mild‐to‐moderate AD. An explanation for this discrepancy might be that we analysed α‐diversities at the class and genus rather than at the species or the amplicon sequence variants (ASVs) level because the V1‐3 region of the 16S rRNA gene operon that we used for analyses does not allow proper resolution at the species level or beyond.

Analysis of the β‐diversity showed distinct clusters for severe AD patients for both fungal and bacterial communities. We thus confirm previous observations indicating an association between bacterial communities and severity of AD. However, in contrast to previous investigations, the microbiota in lesional skin was not different from the one of unaffected skin. , Last, we did not observe an alteration of the skin microbiota within 4 weeks for both AD patients and HC. A previous study analysing skin microbiota in different disease states has not seen any fluctuations in the mycobiota but observed an accumulation of S. aureus during disease flares.

Our study has limitations. First, due to the relatively small study population we were not able to analyse covariates such as treatment, skin care regimens or type of the disease (e.g. head and neck AD). The effect of topical corticosteroids could also not be analysed because all AD patients had been recently treated with them. For dupilumab, an anti‐IL‐4/IL‐13 antibody, 5/7 patients with severe AD and only 1/9 patient with mild‐to‐moderate AD used this drug. Thus, when comparing the two groups, we basically compared severe AD to mild‐to‐moderate AD and the results were mostly corresponding (data not shown). Further, all participants were allowed to use their regular skin care habits. Cosmetics might have an influence on the skin microbiome, which was not controlled and investigated in this study. However, the sample size was big enough to analyse the different severities of the disease with statistically confident results. Second, we endeavoured to use a cohort that reflects the natural course of the disease, which resulted in heterogeneous data with some patients having lesions at one skin site, whereas others did not, making it difficult to compare to other studies. Third, to be able to draw a conclusion on the impact of time in AD, either more time points would be needed or a classification into baseline, flare and postflare should be done. Fourth, with the focus on Malassezia spp., we sampled typical head and neck AD sites. It would be interesting for future studies to include more body sites, for example from the lower part of the body.

In summary, we have observed that (i) non‐Malassezia fungi were increased in patients with severe AD resulting in a higher diversity of the mycobiota, (ii) staphylococci dominated the bacterial communities in AD and (iii) fungal and bacterial skin microbiota of severe AD patients significantly varied from mild‐to‐moderate AD patients and HC. Thus, we conclude that the skin communities of HC and patients with mild‐to‐moderate AD are overall similar with some distinctions in the bacterial communities, that is more staphylococci and less cutibacteria in mild‐to‐moderate AD. However, severe AD is associated with a pronounced dysbiosis of the microbiota with an increase of potential infectious agents. Larger studies are needed to investigate covariates such as treatment with dupilumab or skin care regimens.

Appendix S1. Supporting Information.

Click here for additional data file.

Figure S1 Relative abundance of fungal classes and Malassezia species for each sample.

Click here for additional data file.

Figure S2. Effect size plots of the most abundant fungal genera.

Click here for additional data file.

Figure S3. Relative abundance of bacterial classes and of Staphylococcus species for each sample.

Click here for additional data file.

Figure S4. Effect size plots of the most abundant bacterial genera.

Click here for additional data file.

Figure S5. Spearman correlations of Staphylococcus aureus and SCORAD.

Click here for additional data file.

Figure S6. Fungal ß‐diversities of lesional and unaffected skin in AD patients displayed in principal coordinates analysis (PCoA) plots.

Click here for additional data file.

Figure S7. Bacterial ß‐diversities of lesional and unaffected skin in AD patients displayed in principal coordinates analysis (PCoA) plots.

Click here for additional data file.

Figure S8. Fungal ß‐diversities at each visit displayed in principal coordinates analysis (PCoA) plots.

Click here for additional data file.

Figure S9. Bacterial ß‐diversities at each visit displayed in principal coordinates analysis (PCoA) plots.

Click here for additional data file.