https://orcid.org/0000-0002-2361-7228
Figures
Abstract
Background
Seborrheic dermatitis (SD) is a chronic inflammatory skin disease with a multifactorial aetiology. Malassezia yeasts have been associated with the disease but the role of bacterial composition in SD has not been thoroughly investigated.
Objectives
To profile the bacterial microbiome of SD patients and compare this with the microbiome of individuals with no inflammatory skin disease (controls).
Methods
This was a cross sectional study embedded in a population-based study. Skin swabs were taken from naso-labial fold from patients with seborrheic dermatitis (lesional skin: n = 22; non-lesional skin SD: n = 75) and controls (n = 465). Sample collection began in 2016 at the research facility and is still ongoing. Shannon and Chao1 α- diversity metrics were calculated per group. Associations between the microbiome composition of cases and controls was calculated using multivariate statistics (permANOVA) and univariate statistics.
Results
We found an increased α-diversity between SD lesional cases versus controls (Shannon diversity: Kruskal-Wallis rank sum: Chi-squared: 19.06; global p-value = 7.7x10-5). Multivariate statistical analysis showed significant associations in microbiome composition when comparing lesional SD skin to controls (p-value = 0.03;R2 = 0.1%). Seven out of 13 amplicon sequence variants (ASVs) that were significantly different between controls and lesional cases were members of the genus Staphylococcus, most of which showed increased composition in lesional cases, and were closely related to S. capitis S. caprae and S. epidermidis.
Citation: Sanders MGH, Nijsten T, Verlouw J, Kraaij R, Pardo LM (2021) Composition of cutaneous bacterial microbiome in seborrheic dermatitis patients: A cross-sectional study. PLoS ONE 16(5): e0251136. https://doi.org/10.1371/journal.pone.0251136
Editor: Thomas L. Dawson, Skin Research Institute Singapore, SINGAPORE
Received: December 23, 2020; Accepted: April 20, 2021; Published: May 24, 2021
Copyright: © 2021 Sanders et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: Data cannot be shared publicly as the data involves sensitive research participants data of The Rotterdam Study. Interested researchers may contact our data management team (secretariat.epi@erasmusmc.nl) for access to sensitive data.
Funding: The Rotterdam Study is funded by Erasmus Medical Center and Erasmus University Rotterdam; Netherlands Organization for the Health Research and Development (ZonMw); the Research Institute for Diseases in the Elderly (RIDE); the Ministry of Education, Culture and Science, the Ministry for Health, Welfare and Sports, the European Commission (DG XII), and the Municipality of Rotterdam. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
Abbreviations: SD, Seborrheic dermatitis; OTUs, operational taxonomic units; PCA, principal component analysis; ASVs, amplicon sequencing variants
Introduction
Seborrheic dermatitis (SD) is a common chronic relapsing inflammatory skin disease occurring mainly on the sebum-rich areas of the face, scalp and chest. The disease spectrum is heterogeneous ranging from a mild form of scaling confined to the scalp to severe erythema and scaling of the scalp, face and trunk. The aetiology of SD is complex. Host factors such as skin barrier dysfunction, genetic susceptibility and immune status may all influence disease risk [1–3]. Also, several life style and environmental factors increase risk for SD. In recent epidemiological studies this has been shown for male gender, white skin colour, winter climate, generalized dry skin and specific diets [4, 5].
Since the condition was first described in 1887 [6], it is often assumed that overgrowth of species of the genus Malassezia cause and maintain SD. The Malassezia yeast are more abundant on affected skin, and treatments targeting this yeast (e.g. ketoconazole) can effectively treat SD [7]. However, the genus Malassezia is also part of the healthy microbiome and similar species of Malassezia in SD patients and controls have been found. There might be a causal relationship between Malassezia and SD, but it is not the only factor [8]. A few more recent studies have shown that next to a higher abundances of Malassezia yeast, SD patient might also have bacterial dysbiosis on the scalp [9, 10]. One other study investigated the microbiome on the face of 24 patients with SD and found a higher density of Acinetobacter, Staphylococcus and Streptococcus in lesional skin when compared with non-lesion skin [11]. As this study did not include healthy controls and these bacteria are also part of the normal skin microbiome, it is not clear whether they are indeed associated with SD. Also, members of the Staphylococcus genus can act as pathogens and thus it would be of interest to investigate which specific species are relevant in SD. In addition, as for any observational studies, initial findings need to be replicated in an independent cohort.
As stated above, most studies have focused on the role of Malassezia or were based on small series of cases with no controls, which are necessary to control for ‘normal’ variation in the microbiome. Here we profiled the bacterial microbiome of the nasolabial fold of participants with SD in a middle aged and elderly population based study in the Netherlands (Rotterdam Study) and compared this to the same skin areas of individuals without SD. By using exact amplicon sequence variants (ASVs) instead of operational taxonomic units (OTUs), we were able to zoom in at species level and made our results easier to compare [12].
Materials and methods
Study participants and sample collection
This is a cross-sectional study conducted within the Rotterdam Study. The Rotterdam Study started in 1990 and is an ongoing prospective population based cohort study in a middle aged and elderly population in the Ommoord district of Rotterdam, The Netherlands. All residents aged 45 years and above were invited to participate. At this point, the Rotterdam Study compromises over 15,000 subjects of predominantly North-European descent. Details of the study have been published before [13]. The Rotterdam Study has been approved by the Medical Ethics Committee of the Erasmus MC (registration number MEC 02.1015) and by the Dutch Ministry of Health, Welfare and Sport (Population Screening Act WBO, license number 1071272-159521-PG). The Rotterdam Study Personal Registration Data collection is filed with the Erasmus MC Data Protection Officer under registration number EMC1712001. The Rotterdam Study has been entered into the Netherlands National Trial Register (NTR; www.trialregister.nl) and into the WHO International Clinical Trials Registry Platform (ICTRP; www.who.int/ictrp/network/primary/en/) under shared catalogue number NTR6831. All participants provided written informed consent to participate in the study and to have their information obtained from treating physicians.
In 2016, collection of skin swabs of the nasolabial fold were introduced and are still ongoing. Participants were not allowed to use any skincare products on the day of the examination. Participants that did use skincare products were excluded (n = 108). Participants who used antibiotics in the past 6 months were also excluded (n = 126). Data collection for this research included participants with swabs until April 2018.
Seborrheic dermatitis was diagnosed by dermatology-trained physician during a full body skin examination. The SD diagnosis was based on a greasy scaling, erythema and a characteristic distribution in areas rich in sebaceous glands. Since the microbial swabs were taken from the nasolabial folds in the Rotterdam Study, participants with SD and involvement of the nasolabial fold were considered lesional cases. Participants with SD, but without involvement of the nasolabial fold were considered non-lesional cases. Participants in the Rotterdam Study without SD or any other skin diseases were considered as controls. No sample size calculation was performed.
Swab collection, DNA extraction and 16S rRNA gene polymerase chain reaction amplification and sequencing
A cotton wool swab, pre-moistened with three drops of 0.9% NaCl fluid, was rubbed forth and back, in parallel position along the nasolabial fold (50 times, during 30 seconds). Samples were stored at -20°C for a maximum period of three hours and then at -140°C before processing. Negative swabs (air swabs) were taken for quality control of the sequencing.
DNA was extracted using the DNA Extraction Kit on the Arrow pipetting instrument (DiaSorin S.P.A., Saluggia, Italy). Swabs were treated with DNA Pre-treatment Buffer 2 and Proteinase K for 30 minutes at 56°C, followed by DNA isolation on the Arrow instrument according to the manufacturer’s protocol in batches of 12 samples per run. DNA concentration was measured using the Quant-iT PicoGreen dsDNA Assay Kit (Thermo Fisher Scientific, Waltham, MA) and DNA was stored at -20°C.
The V1 to V3 variable regions of the 16S rRNA gene were amplified using the 27F-519R primer pair and dual indexing as previously described [14, 15]. Amplicons were normalized and pooled in batches (total number of swabs per sequencing run: 300). The pools were purified using Agencourt AMPure XP (Beckman Coulter Life Science, Indianapolis, IN) and the quantity of the pool was assessed using the Quant-iT PicoGreen dsDNA Assay Kit. PhiX Control v3 library (Illumina Inc., San Diego, CA) was spiked into (~10%) the pool prior to sequencing on an Illumina MiSeq sequencer (MiSeq Reagent Kit v3, 2 x 300 bp). Samples were sequenced in three batches.
Bioinformatic analysis
Raw reads from Illumina MiSeq were demultiplexed using a custom script to separate sample fastq files based on the dual index. Primers, barcodes and heterogeneity spacers were trimmed off using tagcleaner v0.16 [16]. Trimmed fastq files were further processed using the DADA2 pipeline [17]. In contrast to standard clustering of reads that cluster sequence reads based on a similarity threshold (usually 97%), and assign them to taxa to generate operational taxonomic units (OTUs), DADA2 allows to analyse exact amplicon sequence variants (ASVs) as the unit of analysis without imposing dissimilarity thresholds [17]. This is done by inferring the biological sequences before incorporating amplification and sequencing errors. Parameters used to run DADA2 are detailed in the S1 File.
Data analysis
Demographic characteristics of participants included in the study were presented as proportions and compared between cases and controls using chi-squared statistics. We excluded ASVs belonging to contaminant phyla as detailed in Salter et al. [18]. Further, we included ASVs that were present in at least 5% of the total dataset and had a minimum of 1000 reads based on the percentile count distribution. Frequency plots of most abundant phyla in the dataset were calculated using relative abundances. Most abundant genera were plotted using different prevalence thresholds from (5 to 100%) using a heatmap plot.
To compare the profile of bacteria between SD patients against controls, we first calculated α-diversity per sample using Chao1 and Shannon diversity. Chao1 is an estimation of the richness of the sample (aka: how many ASV per group), while Shannon-diversity gives an estimate of the relative distribution of sequences among the ASVs. We tested for differences in Chao1 and Shannon diversity using Kruskal-Wallis tests and adjusted the p-values for pair-wise comparisons. We also calculated Shannon and Chao1 after rarefying all samples to have exactly 10000 counts. These analysis were done to evaluate whether any significant difference between the groups was due to different library sizes (details are presented in the S1 File).
Next, for downstream analysis we transformed the dataset using centered log-ratio approach to account for unequal library sizes. This data transformation is recommended to analyse data derived from Next Generation Sequencing since this data is assumed to be compositional (data conveys only relative information since the total abundance is unknown) [19, 20]. Because the logarithmic normalization does not handle zeros, these needed to be imputed first. The imputation of zeroes can be done under different assumptions (see Palarea-Albaladejo et al. for details [21]). Here, we assumed that any feature (ASV) observed in more than one sample could appear in another sample if sequenced with infinite depth, which can be modelled using Bayesian-Multiplicative replacement of count zeros [21].
To visualize whether the microbiome composition clustered per sample characteristics (here, being a case or control, or other sample characteristics), we used Principal Component Analysis (PCA). Further, we tested for differences in the microbiome composition using permutational multivariate analysis of variance (permANOVA) [22] and adjusting for sex, age (as a continuous variable) and batch variables.
Finally, we performed compositional data analysis to test for differential abundance between specific ASVs for seborrheic (lesional and non-lesional) cases and controls with lesional cases used as the reference groups using the ALDEx2 package [23]. The package uses centered log-ratio transformation of the data to correct for unequal library sizes. We tested for differential composition using a generalized model, while adjusting for sex, age and batch as covariates. Significant p-values were defined with a threshold of p-value = 0.05 and borderline with a p-value of 0.1. Of note, the package considers significant effects when the absolute value of the coefficients are higher than 1. However we considered estimates significant if the p-value was <0.05 after multiple testing (Benjamini-Hochberg Procedure, adjusted p-values are given by the package).
For the analysis of alpha and beta diversities the phyloseq and microbiome packages were used [24] (https://github.com/microbiome/microbiome). The pairwise comparison and graphical display of the alpha-diversity was performed using library ("ggstatsplot") retrieved from https://cran.r-project.org/web/packages/ggstatsplot/index.html.
PermANOVA analysis was done using the function “adonis” from the vegan package (https://cran.r-project.org/web/packages/vegan/index.html). Imputation of zeroes was done using the Z compositions packages. The generation of plots and statistical analysis were done using phyloseq and microbiome packages. All the packages were implemented in Rv.3.5.2 (Eggshell Igloo version). Functions and parameters used for these analysis are detailed in the S1 File.
Phylogenetic analysis
To improve the taxonomic resolution of the ASV that showed differential composition between lesional cases and controls, an exploratory phylogenetic analysis was done. We performed this by aligning the sequence reads of the ASVs that showed differential composition between cases and controls with the partial 16S rRNA sequences from different Staphylococcus species retrieved from GenBank (reference numbers were obtained from Ghebremedhin B. et al.) [25]. After the multiple sequencing alignment, a phylogenetic tree was constructed using the neighbour-joining algorithm. The phylogenetic analysis was implemented in the package MEGA [26]. Details on the programs and parameters for generation of the phylogenetic tree are presented in the S1 File.
Results
Skin swabs were collected from 562 participants. Of these 75 participants were classified as non-lesional cases and 22 as lesional cases. Demographic characteristics are shown in Table 1. Most of the participants were male in both non-lesional and lesional cases. The median age was 53 years (inter-quantile range; IQR: 48–64 years) in the control group, and 56 years (IQR: 49–55 years) for non lesional cases and 68 years (IQR: 49–55 years) for lesional cases. The quantile age distribution was significantly different with lesional cases being older than the controls.
Bacterial community profiles of the naso-labial fold
In total, 15651 amplicon sequencing variants (ASVs) were identified in 562 individuals. The ASVs were assigned to 27 unique phyla (S1 Fig), some of which were clearly contaminant phyla including Planctomycetes and Cyanobacteria. The total number of reads was 12277291 with a median count of 21340 (range: 0–218661) per sample. After quality control that included the filtering of contaminant ASVs and samples with low counts and outliers (S1 File), and including ASVs present in at least 5% of the total dataset, 185 ASVs in 529 samples (controls: 436; non-lesional cases:71; lesional cases:22) remained.
Fig 1 presents the relative frequency of phyla in the filtered dataset averaged per control SD cases groups. As shown in Fig 1, we identified five major phyla in the data set, namely: Actinobacteria, Proteobacteria, Bacteroidetes, Firmicutes and Fusobacteria.
There was an apparent increase in the relative abundance of Firmicutes in lesional cases when compared with controls at phylum level as seen in the plot. However, testing the differential composition of SD lesional, non lesional and controls was not significant. Fig 2 showed that several genera seemed to be increased in the lesional cases. However, after testing, only Corynebacterium_1 was nominally significant (p-value< = 0.05; Aldext text; S2 Table).
The relative abundance of main genera is presented separately, with controls = 0, non-lesional cases = 1 and lesional cases = 2. Only Corynebacterium_1 was nominally significantly different between lesional SD versus controls.
Fig 3 presents a heatmap of ASVs identified at different prevalence estimates in all samples (y-axis) and at different levels of detections per sample; from 0.0001% to 20% (x-axis). As shown in this plot, one member of the genus Cutibacterium (ASV1) was the most prevalent and was found in over 75% of all samples with relative abundances of 3% and higher. The second most prevalent ASV belonged to the genus Staphylococcus (ASV2) with a prevalence of 75% at lower abundances (< = 2%). Other ASVs present in at least 50% of the samples included other ASVs belonging to Staphylococcus species (ASV68, ASV5, ASV18; Staphylococcus epidermidis), Cutibacterium granulosum (ASV7) Peptoniphilus rhinitidis and Anaerococcus.
Fig 3 presents the most prevalent ASVs (y-axis) identified at different detection thresholds (0.001 to 20%). Most common ASVs appear at low abundance per sample (0.01 to 1%).
To investigate whether the participants clustered based on their disease status we used PCA. S2 Fig displays the first two axis of variation in the cutaneous microbiome composition. The first component explained 13% of the variation of the data and the second component explained 7%. However, there was not a clear separation between cases and controls.
Comparison of the microbiome composition between SD and controls
We estimated α-diversity of our data per group (controls; non-lesional and lesional SD) using Shannon-diversity and Chao1 statistics. As shown in Fig 4A, the median Chao1 was significantly higher for lesional SD cases in comparison with both non-lesional cases and controls (Kruskal-Wallis rank sum: Chi-squared: 11.14; global p-value = 0.004; controls vs lesional cases; adjusted p-value = 0.003; lesional cases vs non-lesional cases; adjusted p-value = 0.0012). Shannon α-diversity estimates were also significantly different between lesional cases and controls (Kruskal-Wallis rank sum: Chi-squared: 19.06; global p-value = 7.3x10-5; lesional cases vs controls; adjusted p-value = 2.4x10-4; lesional cases vs non-lesional cases; adjusted p-value = 0.03; non-lesional cases versus controls: adjusted p-value = 0.03). Since alpha diversity parameters are sensitive to differences in sample library size, we evaluated the impact of different filtering thresholds as well as rarefaction (S1 File). As shown in the S3–S5 Figs, alpha diversity estimates were significantly higher in lesional cases when compared with non-lesional cases and controls at different filtering thresholds and after rarefying (S1 File).
Fig 4A presents the violin-plot of Chao1 diversity in controls, lesional and non-lesional cases. Fig 4B presents the violin-plot of Shannon diversity in controls, lesional and non-lesional cases. P-values are adjusted pair-wise comparisons. Overall Kruskal-Wallis test is described in the text.
Multivariate analysis.
We tested for differences in the microbiome composition between SD cases and controls using permANOVA. We found that overall, the composition of the microbiome differed between cases and controls (R2 = 0.1%; p-value = 0.03). In addition, age, sex and batch effects were also significant (age: p-value = 0.009; R2 = 0.6%; sex: p-value: 001; R2 = 2%; batch: p-value = 0.001; R2 = 1%).
Univariate analysis.
We investigated whether there was a difference in the composition of ASVs between SD and controls adjusting for known sex, age and batch variables. We found nominally significant differences in the composition between lesional cases compared against controls for bacteria of the genus Staphylococcus including ASV18 (effect size -3.37; p-value = 0.01), ASV15 (effect size: -2.97; p-value = 0.02); ASV22 (effect size -2.88; p-value = 0.03), ASV13 (effect size -2.64; p-value = 0.04) and ASV5 (effect size -2.41;p-value = 0.05). Other bacteria with differences between lesional and controls included a member of the genus Micrococcus (ASV137), Finegoldia (ASV60), and Lawsonella (ASV32). None of the differences were significant after correcting for multiple testing (Table 2). Cutibacterium (ASV1), which was most prevalent in this cohort, and nominally significant at genus level was only borderline associated with SD (effect size:-0.83, p-value: 0.07).
Phylogenetic analysis of ASV with suggestive differential composition between lesional SD and controls
As shown in the Table 2, seven out of 12 ASVs that were (nominally) significantly different between controls and lesional cases were members of the genus Staphylococcus, most of which showed increased composition in lesional cases in comparison with controls. Of these, none of them could be annotated to species level by the DADA2 pipeline. To better understand the role of this genus in SD, we compared the DNA sequence of the reads with the 16S rRNA sequences available from different strains of Staphylococci available from GenBank (Fig 5). As shown in this figure, the ASVs with differences in composition between lesional cases and controls were closely related to S. capitis (ASV5, ASV13, ASV15; ASV22) S. caprae (ASV) and S. epidermidis (ASV33) and were separate from the from S. aureus sequence.
Fig 5 depicts the phylogenetic relationship between ASVs assigned to Staphylococcus species and Staphylococcus species obtained from GenBank (16S rRNA sequence). This sequences begin with AY*. ASVs closely related to S. capitis belong to the same node (ASV5, ASV13, ASV15 and ASV22). ASV68 appears to belong to S. caprae together with ASV18. ASV63 clustered with S. epidermidis. Sequence of S. aureus appear as outgroup (sequence with most differences).
Discussion
In this large cross-sectional study comparing the cutaneous microbiome of individuals with SD with individuals with no skin disease (controls), we found a significant increase in the α-diversity between lesional cases and controls. Comparing the overall microbiome composition between cases and controls using a multivariate approach also led to significant differences between the groups, although the contribution of skin condition was small in comparison to other sample characteristics such as sex. In the design of the study we assumed that the bacterial microbiome of non-lesional skin would also differ from healthy participants as Soares et al. showed this for people with dandruff [27], but we did not find a clear separation between non-lesional cases and controls. In the univariate analysis we found differences in the composition of members of specific bacteria from the genus Staphylococcus, thought to be members of the normal skin flora, including: S. capitis, S. caprae and S. epidermidis, with most of them being increased in lesional cases compared to controls. No differences were found for S. aureus in this analysis. Other bacteria with an increased abundance in the cases were members of the genus Micrococcus, Finegoldia sp. and Streptococcus sanguinis. The changes found in the univariate analysis were not statistically significant after correcting for multiple testing. Though, especially for the Staphylococcus genus there seems to be a trend, as seen in the comparison between non-lesional and lesional SD groups (Table 2). The principal component analysis did not show a clear separation of microbiome composition between controls and cases. In the previous smaller study by Tanaka et al., a separation between cases and control was suggested, although their graph leaves room for discussion and no statistical test was performed [28]. We replicated the findings that Staphylococcus and Streptococcus have an increase abundance in lesional skin.
In our study we found similar distribution of bacteria from sebum rich areas in our dataset as in previous studies [29]. Using a prevalence of 50% to define a ’core’ microbiome, we found that members of Cutibacterium species (ASV1_ Cutibacterium _NA and Cutinebacterium granulossum) and Staphylococcus species (Staphylococcus epidermidis and Staphylococcus _NA) were the most abundant, which serves as an internal validation of our dataset. Of this ‘core’ bacteria only one member of the Cutibacterium species seemed to have differences in their ‘abundance’ between SD cases and controls as shown in the univariate analysis (Table 2). The other ASVs with suggestive differences in the abundance between cases and controls were not amongst the core set, which suggests that the most differences between lesional seborrheic dermatitis and controls occur for bacteria that have relatively low (relative) abundances between 5% to 10%.
We found differences in the relative composition of members of the Staphylococcus genus, which supports previous findings of a differential composition of this genus in seborrheic dermatitis based on smaller series of cases [11]. Staphylococcus species include different species known to be part of normal skin flora, such as S. epidermidis, S. capitis and S. caprae, as well as pathogens such as S. aureus. Because in most studies on seborrheic dermatitis the analysis of bacterial microbiome is done at genus level, it is difficult to compare our results with these from previous studies since we chose to analyse exact ASV instead of grouping ASVs into genera. In a recent study an increased abundance of S. epidermidis in patients with impaired skin barrier was found [30]. Because commensal bacteria can also act as opportunistic pathogens under special circumstances (e.g.: impaired barrier function), the previous study as well as our results supports the hypothesis that commensal bacteria rather than pathogens may be associated with SD.
The higher alpha diversity in SD cases (both lesional and non lesional) in comparison to the controls is in contrast with studies on bacterial microbiome in other chronic inflammatory skin conditions such as atopic dermatitis, where a decreased alpha diversity in lesional skin has been observed. Although both diseases are associated with a skin barrier impairment, it is likely that the role of microbial dysbiosis is different for the two skin conditions. In atopic dermatitis, alpha diversity is lower due to the predominance of one bacterium, namely: S. aureus, which is considered more of a pathogen than a commensal [31]. Higher S. aureus colonization has been associated with immune dysfunction and a decrease of the production of antimicrobial peptides, which will favour the overgrowth of this bacteria while inhibiting the growth of competing bacteria [32]. We looked at the total number of 47 ASVs related to Staphylococcus species in our data and did a phylogenetic analysis with all sequences available at GenBank for 16S rRNA for a total of 75 sequences. Interestingly, only one of 47 ASVs in our data was phylogenetically closer to S. aureus than to the rest (ASV109: S6 Fig), which suggests that S. aureus is not a major player in SD. It could also be argued that an affected nasolabial skin could be a niche for bacteria from nearby areas such as mouth (e.g. Streptococcus sanguinis) or nose (Peptonilus rinitidis). ASVs associated with these bacteria seemed to be increased in SD (when compared with controls), although the differences were not statistically significant. This also may explain the increase in the diversity seen in our dataset. Thus, one could argue that in SD subjects the increased diversity is secondary to an impaired skin barrier, which leads to the colonization of otherwise commensal bacteria.
As mentioned above, we chose exact ASVs as the unit of analysis instead of OTUs in an attempt to gain taxonomic resolution beyond genus level, since the pipeline can annotate 1 to 2 nucleotides difference between ASVs. However, because the DADA2 bioinformatic pipeline uses a strict threshold to annotate ASVs to bacteria from reference databases such as SILVA (100% match; in other words identical sequences), many of the ASVs could not be annotated to species level, and were left as ‘unknown’ species. In most studies, this issue is solved by grouping ASVs or OTUs into genera as unit of analysis, which may lead to spurious results, because different species with different roles in skin biology or pathology (e.g.: S. aureus versus S. epidermidis) are then assumed to be the same. With the methods we used, we could have agglomerated similar ASVs into species level (e.g. different ASVs associated with S. epidermidis into one single S. epidermidis), but due to the high proportion of ASVs that were unknown at species level, we would not have enough data to analyse. With this, we showed that trying to separate species from the same genus may be valuable to understand the role of specific bacteria in SD. Because whole genome sequencing is still expensive for large studies, especially when analysing microbiome from low to medium biomass, a better annotation of reference datasets and perhaps sequencing of another region or even the whole gene could improve the resolution in future studies.
Our study has some limitations. First, this study includes a middle aged and elderly population, and as the microbiome composition changes with age, this might limit the generalizability to younger patients. Second, the bacterial biomass of the skin is lower than that from other human niches [33]. This has implications for sequencing since samples with low to medium levels of extracted DNA are more prone to contamination by other bacteria, which can lead to spurious results. Here we excluded obvious contaminant phyla based on previous publications and filtered ASVs based on prevalence. Third, although our sample of cases and controls is much larger than most studies on skin microbiome, the number of lesional cases was probably too low in comparison with the number of controls, which had much larger variation and made it difficult to separate cases from controls in the PCA. The low density of skin bacterial DNA also creates problems in the statistical analysis since the data were highly sparse. We used methods that alleviated this issue (imputing zeroes), but still the analysis may have been too conservative to survive multiple testing. Further, most methods assume that rare species are not important in determining variation but this may not be true as shown here in our study. To alleviate this, larger sample sizes will be needed to confirm our results.
To conclude, we found differences in the microbiome composition in a large sample of patients with SD when compared with controls. There was a trend towards higher prevalence of members of the genus Staphylococcus, most of which showed increased composition in lesional cases, and were closely related to S. capitis S. caprae and S. epidermidis. The associations were not significant when correcting for multiple testing. Future studies are required to validate these findings.
Supporting information
S1 Fig. Frequency of 27 unique phyla presented.
The graphs present the taxa prevalence against total counts, with each dot representing a different taxa.
https://doi.org/10.1371/journal.pone.0251136.s001
(TIF)
S2 Fig. Principal component analysis (PCA) of microbiome composition (the centered log- ratios).
No clear separation of microbiome composition was observed between controls (grey dots) and cases separated in non-lesional (blue dots) and lesional SD cases (green dots)
https://doi.org/10.1371/journal.pone.0251136.s002
(TIF)
S3 Fig. Boxplot of α-diversity with minimal filtering.
S3A Fig presents the violin-plot of Chao1 diversity in controls, lesional and non-lesional cases. S3B Fig presents the violin-plot of Shannon diversity in controls, lesional and non-lesional cases. P-values are adjusted pair-wise comparisons.
https://doi.org/10.1371/journal.pone.0251136.s003
(TIF)
S4 Fig. Boxplot of α-diversity before 5% prevalence filtering.
S4A Fig presents the violin-plot distribution of Chao1 diversity in controls, lesional and non-lesional cases. S4B Fig presents the violin-plot of Shannon diversity in controls, lesional and non-lesional cases. P-values are adjusted pair-wise comparisons.
https://doi.org/10.1371/journal.pone.0251136.s004
(TIF)
S5 Fig. Boxplot of α-diversity after rarefying.
S5A Fig presents the violin-plot of Chao1 diversity in controls, lesional and non-lesional cases. S5B Fig presents the violin-plot of Shannon diversity in controls, lesional and non-lesional cases. P-values are adjusted pair-wise comparisons.
https://doi.org/10.1371/journal.pone.0251136.s005
(TIF)
S6 Fig. Phylogenetic tree of the ASVs belonging to the genus Staphylococcus.
The figure depicts the phylogenetic relationship between ASVs assigned to Staphylococcus species and Staphylococcus species obtained from GenBank (16S rRNA sequence). This sequences begin with AY*. Only ASV109 was phylogenetically closer to S. aureus.
https://doi.org/10.1371/journal.pone.0251136.s006
(TIF)
S1 Table. Association analysis between microbiome composition (phylum level) and seborrheic dermatitis (with categories controls, non-lesional cases and lesional cases.
a Effect sizes and p-values of the comparison: controls versus lesional cases. b Effect sizes and p-values of the comparison: non-lesional cases versus lesional cases.
https://doi.org/10.1371/journal.pone.0251136.s007
(DOCX)
S2 Table. Association analysis between microbiome composition (genus level) with seborrheic dermatitis (with categories controls, non-lesional and lesional cases) and with lesional cases as reference groups.
a Effect sizes and p-values of the comparison: controls versus lesional cases. b Effect sizes and p-values of the comparison: non-lesional cases versus lesional cases.
https://doi.org/10.1371/journal.pone.0251136.s008
(DOCX)
References
- 1. Chatzikokkinou P, Sotiropoulos K, Katoulis A, Luzzati R, Trevisan G. Seborrheic dermatitis—an early and common skin manifestation in HIV patients. Acta Dermatovenerol Croat. 2008;16(4):226–30. pmid:19111149
- 2. Lally A, Casabonne D, Newton R, Wojnarowska F. Seborrheic dermatitis among Oxford renal transplant recipients. J Eur Acad Dermatol Venereol. 2010;24(5):561–4. pmid:19888944
- 3. Sanders MGH, Pardo LM, Uitterlinden AG, Smith AM, Ginger RS, Nijsten T. The Genetics of Seborrheic Dermatitis: A Candidate Gene Approach and Pilot Genome-Wide Association Study. J Invest Dermatol. 2018;138(4):991–3. pmid:29203360
- 4. Sanders MGH, Pardo LM, Franco OH, Ginger RS, Nijsten T. Prevalence and determinants of seborrhoeic dermatitis in a middle-aged and elderly population: the Rotterdam Study. Br J Dermatol. 2018;178(1):148–53. pmid:28856679
- 5. Sanders MGH, Pardo LM, Ginger RS, Kiefte-de Jong JC, Nijsten T. Association between diet and seborrheic dermatitis: a cross-sectional study. J Invest Dermatol. 2018. pmid:30130619
- 6. Unna PG. Das seborrhoische Ekzem. Monatsh Prakt Der-matol. 1887;6:827–846.
- 7. Kastarinen H, Oksanen T, Okokon EO, Kiviniemi VV, Airola K, Jyrkka J, et al. Topical anti-inflammatory agents for seborrhoeic dermatitis of the face or scalp. Cochrane Database Syst Rev. 2014;5:CD009446. pmid:24838779
- 8. DeAngelis YM, Gemmer CM, Kaczvinsky JR, Kenneally DC, Schwartz JR, Dawson TL Jr., Three etiologic facets of dandruff and seborrheic dermatitis: Malassezia fungi, sebaceous lipids, and individual sensitivity. J Investig Dermatol Symp Proc. 2005;10(3):295–7. pmid:16382685
- 9. Park T, Kim HJ, Myeong NR, Lee HG, Kwack I, Lee J, et al. Collapse of human scalp microbiome network in dandruff and seborrhoeic dermatitis. Exp Dermatol. 2017. pmid:28094891
- 10. Lin Q, Panchamukhi A, Li P, Shan W, Zhou H, Hou L, et al. Malassezia and Staphylococcus dominate scalp microbiome for seborrheic dermatitis. Bioprocess Biosyst Eng. 2020. pmid:32219537
- 11. Tanaka A, Cho O, Saito C, Saito M, Tsuboi R, Sugita T. Comprehensive pyrosequencing analysis of the bacterial microbiota of the skin of patients with seborrheic dermatitis. Microbiology and Immunology. 2016;60(8):521–6. pmid:27301664
- 12. Callahan BJ, McMurdie PJ, Holmes SP. Exact sequence variants should replace operational taxonomic units in marker-gene data analysis. ISME J. 2017;11(12):2639–43. pmid:28731476
- 13. Ikram MA, Brusselle GGO, Murad SD, van Duijn CM, Franco OH, Goedegebure A, et al. The Rotterdam Study: 2018 update on objectives, design and main results. Eur J Epidemiol. 2017;32(9):807–50. pmid:29064009
- 14. Fadrosh DW, Ma B, Gajer P, Sengamalay N, Ott S, Brotman RM, et al. An improved dual-indexing approach for multiplexed 16S rRNA gene sequencing on the Illumina MiSeq platform. Microbiome. 2014;2(1):6. pmid:24558975
- 15. Meisel JS, Hannigan GD, Tyldsley AS, SanMiguel AJ, Hodkinson BP, Zheng Q, et al. Skin Microbiome Surveys Are Strongly Influenced by Experimental Design. J Invest Dermatol. 2016;136(5):947–56. pmid:26829039
- 16. Schmieder R, Lim YW, Rohwer F, Edwards R. TagCleaner: Identification and removal of tag sequences from genomic and metagenomic datasets. BMC Bioinformatics. 2010;11(1):341. pmid:20573248
- 17. Callahan BJ, McMurdie PJ, Holmes SP. Exact sequence variants should replace operational taxonomic units in marker-gene data analysis. The ISME Journal. 2017;11(12):2639–43. pmid:28731476
- 18. Salter SJ, Cox MJ, Turek EM, Calus ST, Cookson WO, Moffatt MF, et al. Reagent and laboratory contamination can critically impact sequence-based microbiome analyses. BMC Biol. 2014;12:87–. pmid:25387460
- 19. Gloor GB, Wu JR, Pawlowsky-Glahn V, Egozcue JJ. It’s all relative: analyzing microbiome data as compositions. Ann Epidemiol. 2016;26(5):322–9. pmid:27143475
- 20. Gloor GB, Macklaim JM, Pawlowsky-Glahn V, Egozcue JJ. Microbiome Datasets Are Compositional: And This Is Not Optional. Frontiers in Microbiology. 2017;8(2224).
- 21. Palarea-Albaladejo J, Martín-Fernández JA. zCompositions—R package for multivariate imputation of left-censored data under a compositional approach. Chemometrics and Intelligent Laboratory Systems. 2015;143:85–96.
- 22. Kelly BJ, Gross R, Bittinger K, Sherrill-Mix S, Lewis JD, Collman RG, et al. Power and sample-size estimation for microbiome studies using pairwise distances and PERMANOVA. Bioinformatics. 2015;31(15):2461–8. pmid:25819674
- 23. Fernandes AD, Macklaim JM, Linn TG, Reid G, Gloor GB. ANOVA-like differential expression (ALDEx) analysis for mixed population RNA-Seq. PloS one. 2013;8(7):e67019–e. pmid:23843979
- 24. McMurdie PJ, Holmes S. phyloseq: an R package for reproducible interactive analysis and graphics of microbiome census data. PLoS One. 2013;8(4):e61217. pmid:23630581
- 25. Ghebremedhin B, Layer F, König W, König B. Genetic classification and distinguishing of Staphylococcus species based on different partial gap, 16S rRNA, hsp60, rpoB, sodA, and tuf gene sequences. Journal of clinical microbiology. 2008;46(3):1019–25. pmid:18174295
- 26. Tamura K, Stecher G, Peterson D, Filipski A, Kumar S. MEGA6: Molecular Evolutionary Genetics Analysis version 6.0. Mol Biol Evol. 2013;30(12):2725–9. pmid:24132122
- 27. Soares RC, Camargo-Penna PH, de Moraes VC, De Vecchi R, Clavaud C, Breton L, et al. Dysbiotic Bacterial and Fungal Communities Not Restricted to Clinically Affected Skin Sites in Dandruff. Front Cell Infect Microbiol. 2016;6:157. pmid:27909689
- 28. Tanaka A, Cho O, Saito C, Saito M, Tsuboi R, Sugita T. Comprehensive pyrosequencing analysis of the bacterial microbiota of the skin of patients with seborrheic dermatitis. Microbiol Immunol. 2016;60(8):521–6. pmid:27301664
- 29. A Grice E, Kong H, Conlan S, Deming C, Davis J, Young A, et al. Topographical and Temporal Diversity of the Human Skin Microbiome2009. 1190–2 p.
- 30. An Q, Sun M, Qi RQ, Zhang L, Zhai JL, Hong YX, et al. High Staphylococcus epidermidis Colonization and Impaired Permeability Barrier in Facial Seborrheic Dermatitis. Chin Med J (Engl). 2017;130(14):1662–9. pmid:28685715
- 31. Otto M. Staphylococci in the human microbiome: the role of host and interbacterial interactions. Curr Opin Microbiol. 2020;53:71–7. pmid:32298967
- 32. Nakatsuji T, Chen TH, Narala S, Chun KA, Two AM, Yun T, et al. Antimicrobials from human skin commensal bacteria protect against Staphylococcus aureus and are deficient in atopic dermatitis. Sci Transl Med. 2017;9(378). pmid:28228596
- 33. Ferretti P, Farina S, Cristofolini M, Girolomoni G, Tett A, Segata N. Experimental metagenomics and ribosomal profiling of the human skin microbiome. Exp Dermatol. 2017;26(3):211–9. pmid:27623553