Identification of cutaneous fungi and mites in adult atopic dermatitis: analysis by targeted 18S rRNA amplicon sequencing

Atopic dermatitis (AD) patients have an altered skin bacterial community, with an abundance of Staphylococcus aureus associated with flares, highlighting that microbial organisms may be important for disease exacerbation. Despite strong evidence of association between bacterial skin colonisation and AD, very limited knowledge regarding the eukaryotic microbial community, including fungi and ectoparasites, in AD exists. In this study, we compared the skin and nasal eukaryotic microbial community between adult AD patients (n = 55) and non-AD healthy controls (n = 45) using targeted 18S rRNA amplicon sequencing. Analysis was based on the presence or absence of eukaryotic microorganisms. The cutaneous composition of the eukaryotic microbial community and the alpha-diversity differed significantly between AD patients and non-AD individuals, with increased species richness on AD skin. Alpha-diversity and beta-diversity were similar on lesional and non-lesional skin of patients. The ectoparasite Demodex folliculorum and the yeast Geotrichum candidum were significantly more prevalent on the skin of AD patients. The prevalence of D. folliculorum on lesional skin was greater among patients recently treated with topical corticosteroid. Malassezia was one of the most frequently detected genera at all sites, with M. globosa and M. restricta being the most prevalent. M. restricta was under represented in the anterior nares of AD patients as compared to the non-AD control population. Significant differences in the eukaryotic microbial communities were found between AD patients and non-AD individuals, with the most striking finding being the significantly overrepresentation of D. folliculorum on AD skin. Whether D. folliculorum can contribute to skin inflammation in AD needs further investigation.

disease severity [4][5][6][7]. While several studies have examined the bacterial community on skin, our knowledge regarding the eukaryotic microbial community on AD skin is very limited. Two sequencing-based studies have shown that the fungal richness and diversity are greater on AD lesional skin (LS) compared to healthy control skin [8,9]. Malassezia species, especially M. globosa and M. restricta, are the most common and abundant fungal species on skin of AD patients and healthy individuals [8][9][10][11][12]. Malassezia has been implicated in AD pathogenesis, as patients are more often hypersensitive with specific IgE antibodies against Malassezia in comparison to healthy individuals, a subtype of AD called 'head and neck dermatitis' [13][14][15]. Candida is another commensal yeast suggested to contribute to the onset and exacerbation of AD [14], but a higher colonisation rate of AD patients compared with controls has not been found [16]. Furthermore, a wide range of environmental molds, such as Aspergillus and Cladosporum species, can induce type I allergic responses, where greater incidence rates have been reported among atopic patients compared to the general population [17].
In this explorative study, we aimed to compare the presence of skin colonising fungi and parasites on adult AD patients and non-AD healthy individuals, using 18S rRNA gene amplicon sequencing. Furthermore, we aimed to investigate whether topical corticosteroid (TCS) treatment, AD disease severity and AD risk factors (filaggrin gene (FLG) mutations and S. aureus colonisation) were associated with changes in the eukaryotic microbial community.

Results
Samples were collected from skin and the nares of 58 adult AD patients and 46 non-AD healthy individuals, including both lesional (LS) and non-lesional (NLS) skin areas from the AD patients. Three patients were excluded prior to analysis due to low quality of samples, and 55 AD patients were thus included in the final analysis. Demographic and clinical descriptions of the study population are given in Table 1.
Compositional analysis of the eukaryotic microbial communities on skin and in nares was based on the presence or absence of specific organisms, where species presence was defined as ≥10 classified sequence reads within a sample (See the Methods Section for more details).
AD skin samples were collected from distinct anatomical sites depending on the presence of eczema, with LS samples primarily collected at dry skin areas (n = 23) such as the volar forearm and dorsal hand, and NLS samples were primarily collected from moist skin areas of the antecubital crease (n = 45) ( Table 1). All skin samples from healthy individuals were also collected from the antecubital crease. Initial analysis indicated that differences in the microenvironment of the sampling area (i.e. dry, moist or sebaceous skin) did not significantly influence the overall eukaryotic microbial composition on either LS or NLS (Additional file 1: Fig. S2) and differences in anatomical sites being sampled between individuals were therefore considered not to significantly influence subsequent analyses.

Eukaryotic microbial community composition and richness
The eukaryotic microbial community composition on AD skin was significantly different from the composition on healthy control skin, for both AD LS (R = 0.12, adjusted p = 0.001; ANOSIM) and AD NLS (R = 0.10; adjusted p = 0.001; ANOSIM), though the close clustering of some AD and control skin samples indicates some degree of similarity between sites (Fig. 1a). No differences were observed between AD LS and NLS (R = 0.00; ANOSIM).
The eukaryotic microbial richness, defined as the total number of observed species in a sample, was overall higher on skin compared to the anterior nares (Fig. 2). Richness was greater on both AD LS (median species count: 29.0 (IQR: 20.5-41.5)) and AD NLS (31.0 (23.5-43.0)) compared to healthy control skin (22.0 (13.0-34.0)), though after correcting for multiple testing only statistically significant between AD NLS and healthy control skin (AD LS vs control skin: p = 0.03, adjusted p = 0.3; AD NLS vs. control skin: p = 0.002, adjusted p = 0.02, Mann-Whitney U test). No significant difference in richness was observed between AD LS and NLS.
Next, we examined whether TCS treatment within one month prior to sampling could influence the diversity on AD LS, and there we observed a significant variation in community composition between the treated (n = 42) and non-treated patients (n = 13) (R = 0.20, p = 0.006, adjusted p = 0.07; ANOSIM) (Additional file 1: Fig. S3). Also, the species richness on AD LS was greater among patients who had been treated with TCS (30.0 (22.5-45.8)) compared to non-treated patients (20.0 (14.0-32.0)), but the difference was not statistically significant after correcting for multiple testing (p = 0.046, adjusted p = 0.4, Mann-Whitney U test) (Additional file 1: Fig. S3).

Single species overrepresentation on atopic dermatitis skin
Two species, the skin mite Demodex folliculorum and the yeast Geotrichum candidum, were significantly more common on AD skin compared to healthy control skin (

Distributions of Malassezia and Candida
The yeast Malassezia was among the most frequently detected genera at all sample sites (Additional file 1: Table S1), with M. globosa and M. restricta being the most prevalent (Fig. 3). In the anterior nares, M. restricta was significantly more prevalent among healthy individuals (98%) compared to AD patients (67%) (adjusted p = 0.03, Fisher's exact test). Of notice, M. furfur was not detected in any samples.
The yeast Candida was more frequently identified on AD skin (LS: 64% and NLS: 65%) compared to healthy control skin (35%); however, the difference was not statistically significant. We tested if treatment with antibiotics within the last 3 months prior to sample collection was associated with the greater prevalence of Candida on AD skin but found no significant difference in the prevalence of Candida on either LS or NLS between patients, treated or not treated with antibiotics (topical and/or systemic treatment) in the past 3 months (Additional file 1: Table S2).

Eukaryotic microbial communities in relation to clinical aspects of atopic dermatitis
There were no significant differences in either species richness or the general community composition on LS and NLS between patients with mild (SCORAD< 25, n = 18), moderate (SCORAD 25-50, n = 31), and severe (SCORAD> 50, n = 6) AD (Additional file 1: Fig. S6). In addition, there  a: Species presence was defined as 10 ≥ classified reads in a sample. b: Pairwise differences between sample sites were tested using Fisher's exact test and odds ratios were calculated using healthy control skin from non-AD individuals as the reference. c: P-values were corrected for mass-significance using the Bonferroni method. Adjusted p-values < 0.05 were considered as statistically significant. Abbreviations: AD atopic dermatitis, LS lesional skin, NLS non-lesional skin, OR odds ratio, CI confidence interval was no difference in either species richness or overall community composition on LS and NLS between patients with (n = 19) or without (n = 31) FLG mutations (Additional file 1: Fig. S7). A hallmark of AD disease is S. aureus skin colonisation. Inter-species interactions and competition could very well influence S. aureus growth and colonisation; however, no single eukaryotic microbial species was over-or underrepresented with respect to S. aureus colonisation in either LS, NLS or nares. Reduced bacterial diversity, measured by Shannon-index, was not associated with changes in eukaryotic microbial richness on AD skin or in the anterior nares (Additional file 1: Fig. S8).

Discussion
Knowledge of the eukaryotic microbial community in AD is limited [8,9]. however, cutaneous fungi and parasites might be of importance in AD pathogenesis, as seen for the bacterial community [18]. In the present study, analysis of similarity implied that the overall eukaryotic microbial community composition on skin and in nares was significantly different between AD patients and non-AD healthy individuals. The community composition was similar on AD LS and NLS, despite clinical differences between the two sample sites. No significant association was found with either AD disease severity (SCORAD) or carriage of loss-of-function FLG mutations, as previously reported for the bacterial community on AD skin [5,19]. The present study confirms that the eukaryotic microbial richness is significantly greater on AD skin compared to healthy control skin [8,9]. This might be due to regular use of TCS among patients, as TCS treatment within 1 month prior to sample collection was associated with an increased richness on AD LS. Another possible explanation might be that antibacterial therapy alters the microbiome by reduction of bacteria thereby contributing to proliferation and changed virulence characteristics of the remaining microbes such as fungi, similar to what have previously been described for Candida vaginitis and Malassezia folliculitis after antibiotic treatment [20][21][22].
To some extent, the observed variation between AD and healthy control skin might be explained by a significantly higher frequency of the ectoparasite D. folliculorum and the fungus G. candidum on AD skin. D. folliculorum is a common skin mite that predominantly lives in sebaceous areas of the face, where it consumes sebum [23]. Though D. folliculorum is a commensal inhabitant of the skin, increased mite densities and penetration into the dermis have previously been associated with inflammation and skin barrier disruption [24,25]. Also, D. folliculorum has been associated with the inflammatory skin disease rosacea [23,26]. D. folliculorum was identified at a significantly higher prevalence on AD skin (LS as well as NLS) in this study, despite the fact that our samples primarily were collected from moist and dry skin areas. A possible explanation for this could be that AD patients regularly use moisturising lotions as part of their treatment, which might create an advantageous environment for D. folliculorum. This hypothesis can be supported by our observation that D. folliculorum was more prevalent on skin among patients reporting use of TCS, which is in alignment with findings from a study of patients with perioral dermatitis [27]. A second hypothesis could be that the skin conditions in AD itself favour D. folliculorum colonisation, as slightly higher colonisation frequencies have been observed among people with higher skin pH and lower skin hydration [28], as present in AD.
G. candidum, the other species found to be overrepresented on AD skin, is considered a common skin colonising yeast [29,30], However, other DNA sequencing based studies on skin fungal communities, including communities on AD skin, have not reported the presence of this organism [8][9][10][11], which might be due to previously small study populations. Whether the distorted skin ecology in AD, including increased skin pH, M. globosa and M. restricta were some of the most frequently observed fungal spp. at all studied sample sites, which is in accordance with previous findings [8,9], suggesting that these species can be considered as commensals of these anatomical areas. Malassezia is a lipophilic yeast, formerly thought only to reside on seborrheic skin areas (scalp, face and thorax), but after the introduction of molecular based detection methods, it has been isolated from most body sites [10]. There was no significant difference in the prevalence of Malassezia spp. on skin between AD patients and healthy controls, which is somewhat surprising as AD patients have a dry skin, and the genus Malassezia is lipid dependant. A significantly higher prevalence of M. restricta was found in the nares of healthy controls compared to the nares of AD patients. The explanation for this observed reduction of M. restricta in AD nares is not clear and needs further examination.
Candida is a commensal of mucosal membranes that colonises moist skin areas, and it is thus surprising that this genus tended to be more frequently detected in AD patients than healthy individuals, as patients with AD have drier skin. High prevalence of Candida with the potential of causing infections has previously been associated with long-term use or repeated use of antibiotics [31], however, antibiotic treatment was not associated with significant changes in Candida presence on AD skin in the present study.
The major strength of this study is the large cohort, which is by far the largest published to date. Furthermore, AD-diagnosis was verified by a specialised dermatologist, and all patients had active disease. This ensures high quality and clinical relevance of the presented data. A well-known limitation in the field of microbiome research is the difficulties in discriminating between colonising organisms and environmentally derived organisms. Thus, some of the identified species could be environmentally derived contaminants of the human epidermis rather than true colonisers of skin and nares, e.g. Cladosporium and Aspergillus are common sporeforming molds in both indoor and outdoor environments [17]. However, increased sensitivity in AD patients to these molds [17], together with their high prevalence found on AD and healthy control skin [8,10], support a possible role for mold sensitivity and AD inflammation.

Conclusion
Significant differences in the composition and richness of the eukaryotic microbial community were observed between AD patients and non-AD healthy individuals, whereas AD disease severity or filaggrin gene mutations were not associated with changes in the eukaryotic microbial community. An unexpected finding was the increased prevalence of the skin mite D. folliculorum on AD skin, which was associated with recent TCS treatment. Colonisation with Demodex spp. has been associated with folliculitis and rosacea [23,25], and it would thus be clinically relevant to examine if D. folliculorum can contribute to skin inflammation in AD.

Study population
Adult AD patients (n = 58) from the outpatient clinic of the Department of Dermatology, Bispebjerg Hospital (Denmark), were included in the study in January-June 2015. Inclusion criteria were age ≥ 18 years and presence of AD according to U.K. criteria [32] as assessed by a specialised dermatologist. Exclusion criteria were pregnancy, breastfeeding and UV therapy within the last 2 months. Patients were invited to participate in the study regardless of topical of systemic treatment of AD; however, information regarding medical treatment 3 months prior to the sampling timepoint was registered (Table 1). Disease severity was assessed using the Severity Scoring of Atopic Dermatitis (SCORAD) index [33]. Blood samples were obtained for identification of the most common filaggrin gene (FLG) mutations among Caucasians (R501X, 2282del4, and R2447X) [34].
Non-AD healthy individuals (n = 46; > 18 years) were recruited among employees at Bispebjerg Hospital and Statens Serum Institut (Denmark) during the months of February and September 2016. Exclusion criteria were current or previous AD, and daily work in microbiology laboratories. The study population has previously been characterized with respect to the bacterial community composition on skin [5,35]. All methods were carried out in accordance with relevant guidelines and regulations.

Sample collection
Swabs (eSwabs, Copan, Italy) were taken from skin (AD LS, AD NLS, and healthy controls) and from anterior nares. Samples from AD NLS and healthy control skin were taken from the antecubital crease, except for 10 patients who had visible eczema at the site. In these cases, NLS samples were primarily taken from the volar forearm. AD LS samples were collected depending on the location of eczema, but primarily from the volar forearm and dorsal hand. Samples were stored at − 80°C.

DNA extraction and sequencing
DNA was extracted from samples using a MagNa Lyser instrument (Roche, Mannheim, Germany) and the FastDNA® SPIN Kit for Soil (MP Biomedicals, Santa Ana, CA, USA). The V3-V5 region of the 18S rRNA gene was amplified in a two-step PCR using three customised primer sets [36]. Three amplicon regions were used in order to increase sequence diversity and taxonomic resolution during classification. Sequencing was carried out on a MiSeq instrument (Illumina Inc., San Diego, CA, USA) using the v2 reagent kit.

Bioinformatics and statistics
BION v.17.10 (http://box.com/bion) was used for sequence mapping (settings described in Additional file 2: Supplemental methods) and the SILVA SSU database v.128 [37] used for taxonomic classification. Malassezia spp. were classified using a custom-made curated database and the DADA2 pipeline [38] (Additional file 2: Supplemental methods).
Analyses were performed in R v.3.5.1 (The R Foundation for Statistical Computing, Vienna, Austria) using the packages phyloseq v.1.24.2 [39], vegan v.2.5-3 [40], and ggplot2 v.3.1.0 [41]. Taxa tables from the three amplicon sequence sets were merged, using the highest observed read count for each species. Sequences classified as phyla belonging to the plant kingdom or the subphylum Vertebrata were removed. Samples < 5000 reads (Additional file 1: Fig. S1) and samples from patients with incomplete sample sets were excluded.
Differences in eukaryotic microbial communities between groups were evaluated by presence/absence analysis, where species presence was defined as ≥10 reads. First, differences in overall community composition between groups were investigated using Jaccard distances and visualized with principal coordinates analysis (PCoA) plots. The degree of similarity was examined using analysis of similarity test (ANOSIM) [42] with 9999 permutations. Second, differences in species richness between groups were examined using Wilcoxon signed rank test (paired samples), Mann-Whitney U test (unpaired samples), or Kruskal Wallis test (>two unpaired groups). Species richness was compared to bacterial Shannon diversity [5] on AD skin and anterior nares using Spearman's rank correlation test. Third, differences in the frequencies of present species between sample sites were examined using Fisher's exact test (species present in 5-95% of samples were included in the analysis). Fisher's exact test was also used to compare presence of pre-selected spp. between AD treatment groups. Results from each analysis were corrected for multiple testing using the Bonferroni method, and adjusted p-values < 0.05 were considered statistically significant (Additional file 2: Supplemental methods).