The genus Staphylococcus is a known component of atopic dermatitis. In this issue, Byrd et al. used shotgun metagenomic sequencing to analyze the species and strains present at baseline and during flares in pediatric atopic dermatitis. They observed that patients with more mild disease had more Staphylococcus epidermidis detected in flares and that those with severe disease were colonized by dominant clonal Staphylococcus aureus strains. S. aureus strains from patients were applied to mouse skin, and strains from severe flares induced T cell expansion and epidermal thickening. These results highlight that not all Staphylococcus can be considered the same, even at the strain level, when it comes to atopic dermatitis pathogenesis.

The heterogeneous course, severity, and treatment responses among patients with atopic dermatitis (AD; eczema) highlight the complexity of this multifactorial disease. Prior studies have used traditional typing methods on cultivated isolates or sequenced a bacterial marker gene to study the skin microbial communities of AD patients. Shotgun metagenomic sequence analysis provides much greater resolution, elucidating multiple levels of microbial community assembly ranging from kingdom to species and strain-level diversification. We analyzed microbial temporal dynamics from a cohort of pediatric AD patients sampled throughout the disease course. Species-level investigation of AD flares showed greater Staphylococcus aureus predominance in patients with more severe disease and Staphylococcus epidermidis predominance in patients with less severe disease. At the strain level, metagenomic sequencing analyses demonstrated clonal S. aureus strains in more severe patients and heterogeneous S. epidermidis strain communities in all patients. To investigate strain-level biological effects of S. aureus, we topically colonized mice with human strains isolated from AD patients and controls. This cutaneous colonization model demonstrated S. aureus strain–specific differences in eliciting skin inflammation and immune signatures characteristic of AD patients. Specifically, S. aureus isolates from AD patients with more severe flares induced epidermal thickening and expansion of cutaneous T helper 2 (T H 2) and T H 17 cells. Integrating high-resolution sequencing, culturing, and animal models demonstrated how functional differences of staphylococcal strains may contribute to the complexity of AD disease.

With an increasing appreciation of functional differences between strains within a single species, we performed shotgun metagenomic sequencing of AD patient skin samples to capture the full genetic potential and strain-level differences of the skin microbiome throughout the course of the disease. We confirmed an increase of staphylococcal species during disease flares in our cohort and more deeply explored the S. aureus and Staphylococcus epidermidis strain diversity of each patient. To test the functional consequence of strain-level differences between patients, we isolated staphylococcal strains from patients and healthy controls and investigated the cutaneous and immunologic effects when applied topically in a mouse model.

In addition to host genetic susceptibility, the relationship between AD and skin bacteria is well recognized clinically. Patients with AD are often treated with varying combinations of antimicrobial approaches (for example, antibiotics and dilute bleach baths) and anti-inflammatory or immunosuppressive medications ( 5 ). The efficacy of these antimicrobial treatments is associated with decreases in staphylococcal relative abundances ( 6 , 7 ). Staphylococcus aureus commonly colonizes AD skin and has been studied using colony-counting, sequencing typing methods [for example, pulsed-field gel electrophoresis and SpA (S. aureus protein A) typing] of selected cultivated isolates ( 8 – 11 ), and more recently, amplicon-based (marker gene) sequencing of the 16S ribosomal RNA gene ( 6 , 7 , 12 , 13 ). However, sequence typing and amplicon sequencing methods are unable to distinguish between genetically distinct strains, as determined by whole-genome sequencing ( 14 , 15 ). By contrast, shotgun metagenomic sequencing of skin samples from healthy individuals provided deeper resolution and demonstrated the multiphyletic composition of commensal Staphylococcus ( 16 ).

Atopic dermatitis (AD; eczema) is a common inflammatory skin disorder in industrialized countries, affecting 10 to 30% of children ( 1 ). Patients with AD suffer from chronic, relapsing, intensely itchy, and inflamed skin lesions and have an increased likelihood of developing asthma and/or hay fever ( 2 ). AD is a complex multifactorial disease in which epidermal barrier impairment, type 2 immunity, and skin microbes are each thought to potentially play a causative role ( 1 ). More than 30 susceptibility loci have been associated with AD, including mutations in the gene encoding the skin barrier protein filaggrin (FLG) ( 3 ) and genes linked to the immune system ( 4 ).

RESULTS

Bacterial communities shift during AD disease progression To examine the relationship between the skin microbiota and AD, 11 children with moderate to severe AD and 7 healthy children were recruited to the National Institutes of Health Clinical Center between June 2012 and March 2015 (tables S1 and S2). Because AD has a chronic relapsing course, patients were sampled at stable/typical disease state (baseline), worsening of disease (flare), and 10 to 14 days after initiation of treatment using a combination of skin-directed therapies (post-flare). Because the use of topical medications on AD skin alters skin bacterial communities (6, 7), baseline samples were defined as those collected from subjects in their routine disease state who refrained from using skin-directed antimicrobial and anti-inflammatory treatments for 7 days, a duration of time determined on the basis of prior findings (6). Flares were defined as time points when patients experienced worsening in the clinical severity of their typical AD, had not used skin-directed antimicrobial and anti-inflammatory treatments for 7 days, and did not have clinical skin infection (for example, yellow crusts or pustules). At each time point, disease severity was determined with objective SCORAD (SCORing Atopic Dermatitis), a validated clinical severity assessment tool (17–19). Subjects were sampled bilaterally at sites of disease predilection: the inner elbow [antecubital crease (Ac)] and behind the knees [popliteal crease (Pc)], along with five additional sites to investigate defined areas with different skin physiologies (fig. S1). Because of the clinical severity of their AD, 6 of the 11 patients experienced exacerbations of their skin disease with the 7-day skin preparation regimen and could not provide baseline time point samples, reflecting the spectrum of the natural course of AD. Because the skin microbial dysbiosis during AD flares was of greatest interest, most of the analyses focused on comparisons between flare and post-flare time points. In total, we performed shotgun metagenomic sequencing of 422 samples, generating 191 giga–base pairs (Gbp) of microbial sequence data from 27 AD patient visits and 7 healthy control visits (table S3). During patient flares, AD disease severity was significantly elevated as indicated by higher mean objective SCORAD (38 ± 2.9) as compared to baseline (9.4 ± 1.6; P < 4.5 × 10−4) and post-flare (11 ± 1.6; P < 2.8 × 10−6) (Fig. 1A). Fig. 1. Bacterial communities shift during AD disease progression. (A) Objective SCORAD for patients at baseline (n = 5), flare, and post-flare (n = 11). Higher SCORAD corresponds to more severe disease. ***P < 0.001, with nonparametric Wilcoxon rank-sum test. (B) Mean Shannon diversity ± SEM in controls and AD disease states. Colors correspond to disease state. Vf, volar forearm; Ic, inguinal crease; Fh, forehead; Oc, occiput; Ra, retroauricular crease. (C) Shannon diversity versus objective SCORAD for AcPc of AD patients. Pearson partial correlation (adjusting for disease state). (D) Mean relative abundance of bacterial genera in AcPc for controls and AD disease states. (E) Mean relative abundance of predominant genera in AcPc for disease states. Statistical significance based on paired Wilcoxon test and Bonferroni correction. F, flare; PF, post-flare. (F) Proportion of Staphylococcus versus objective SCORAD for AcPc of AD patients. Pearson partial correlation (adjusting for disease state). To compare the microbial community composition across time points, we mapped microbial reads to a multikingdom reference database. As seen in healthy adults (16, 20, 21), Bacteria was the most predominant kingdom across time points and body sites (fig. S2 and table S4). Malassezia species, particularly Malassezia restricta and Malassezia globosa, predominated the fungal communities (fig. S3 and table S5), and eukaryotic DNA viral communities were mostly polyomaviruses or papillomaviruses depending on the individual (fig. S4). No significant differences in the fungal or viral components over time were identified; therefore, we focused on bacterial communities that demonstrated the greatest shifts in this cohort (fig. S5 and table S6). We first determined the Shannon diversity index, an ecological measure of richness (total number of bacterial species) and evenness (relative proportion of the bacterial species), to evaluate the overall community structure/composition across body sites and time points. During flares, Ac and Pc, which are the sites of AD predilection, exhibited a marked reduction in Shannon diversity compared to baseline, post-flare, and healthy controls, a trend observed to a lesser extent across other sites (Fig. 1B). Because changes in bacterial diversity were most pronounced at sites of disease predilection and Ac and Pc have similar microbial communities (21), we averaged these sites per subject and used the composite “AcPc” for subsequent analyses. Similar to our previous analysis of microbial diversity in an AD patient cohort (6), the partial correlation between objective SCORAD and Shannon diversity, adjusting for disease state, was significantly inversely correlated (r = −0.58, P = 4.5 × 10−4) (Fig. 1C), indicating that reduced skin bacterial diversity corresponds to worse disease severity, primarily at sites of disease predilection (fig. S5A). To determine which taxa were contributing to the loss of diversity, we compared the relative abundances of the most prominent taxa (Fig. 1D and fig. S5B). Of the four most prominent genera in the AcPc, only Staphylococcus was significantly increased in flares (45 ± 10.2%) as compared to post-flares (9.2 ± 2.4%; P < 0.0078) and healthy controls (6.6 ± 4.1%; P < 0.033) (Fig. 1E). This increase in Staphylococcus relative abundances was positively correlated with objective SCORAD (r = 0.67, P < 8.1 × 10−6) (Fig. 1F), indicating that severe AD was associated with higher staphylococcal relative abundances at sites of disease predilection. In addition, there was a positive correlation for the forehead, retroauricular crease, and volar forearm (fig. S5C), sites that can be affected in more severe disease. However, differences in Corynebacterium, Propionibacterium, and Streptococcus relative abundances between flares and post-flares were not statistically significant (Fig. 1E).

AD flare severity is linked with specific staphylococcal species To further examine the positive correlation between Staphylococcus and AD disease course observed in this study and in prior studies (22), we identified the relative abundances of staphylococcal species including S. aureus, S. epidermidis, Staphylococcus hominis, and Staphylococcus capitis (Fig. 2A and fig. S6). Only relative abundances of S. aureus were significantly increased from flares (28 ± 8.8%) to post-flares (2.3 ± 0.8%; P < 0.014) (Fig. 2B). Although S. epidermidis relative abundances were also higher during flares (13 ± 5.4%) as compared to post-flares (3.7 ± 1.4%), results did not reach statistical significance. For all patients, relative abundances of S. aureus were positively correlated with objective SCORAD (r = 0.73; P < 1.0 × 10−7), whereas S. epidermidis was not correlated (Fig. 2C and fig. S7). This association between S. aureus and AD severity (23) has been observed in prior studies. Neither S. hominis nor S. capitis demonstrated significant shifts in relative abundances between time points (Fig. 2B) or was correlated with disease severity (fig. S7). Fig. 2. Staphylococcal species increase during AD disease flare. (A) Mean relative abundance of staphylococcal species within the total bacterial population in AcPc of AD patients and controls. (B) Mean relative abundance of most abundant Staphylococcus species in AcPc for disease states. Statistical significance based on paired Wilcoxon test. (C) Correlation of S. aureus (left) and S. epidermidis (right) mean relative abundance and objective SCORAD for AcPc of patients. Pearson partial correlation (adjusting for disease state). (D) Comparison of S. aureus to S. epidermidis relative abundance by patient for all sites. The patient’s objective SCORAD is indicated in the parenthesis. Shape corresponds to physiological characteristic of the body site, color to the predominant species, and size to the magnitude of disease severity (objective SCORAD). Patients at the top row have a higher predominance of S. epidermidis, whereas patients at the bottom row are S. aureus–predominant. To further explore the relationship between disease severity and staphylococcal species, we sorted the patients by their objective SCORAD and plotted the relative abundances of S. aureus and S. epidermidis at flare (Fig. 2D). We observed a trend whereby patients with more severe AD flares (objective SCORAD, 45 ± 3.0) had higher relative abundances of S. aureus (Fig. 2D, bottom row; fig. S8; and table S7). In contrast, patients with less severe AD flares, as well as lower objective SCORAD (31 ± 1.9; P < 0.004, in comparison to the more severe flares), had higher relative abundances of S. epidermidis (Fig. 2D, top row) across sampled sites. Specifically, more severe AD flare patients had relative abundances of 34 ± 8.7% S. aureus with 7.4 ± 4.2% S. epidermidis, and less severe AD flare patients had relative abundances of 3.8 ± 1.7% S. aureus with 13 ± 3.9% S. epidermidis averaged across all sites during flare. The range of S. aureus relative abundances based on sequencing was variable: 3 of 11 patients had no S. aureus on their skin, and 3 of 11 patients had relative abundances of S. aureus on their skin exceeding 50%, similar to prior studies of S. aureus relative abundances on AD skin (6, 24, 25). To compare these metagenomic results with more traditional studies, we cultured bacteria from skin and nares swabs collected concurrently with genomic samples. Cultures of S. aureus from skin clinical samples correlated with the microorganism detection by sequencing. Notably, two less severe AD flare patients were culture-positive for S. aureus only in their nares, a common site of carriage. The S. aureus culture-positive rates in this cohort were consistent with those in other studies (6–8, 11–13). The genomic analyses were internally consistent with cultivation results, and both supported the strong association between AD disease severity and S. aureus.

S. aureus strains in AD demonstrate monoclonality Although the differential association of S. aureus and S. epidermidis with AD severity defined an intriguing feature of disease heterogeneity, the underlying strain communities of these species during the disease course remained unknown. Two alternative scenarios could underlie microbial shifts in a disease flare: (i) All strains equally increase in relative abundance, or (ii) a particular strain(s) blooms and drives the increase. This distinction is important because individual strains may exhibit functional differences. In previous studies, this question could not be addressed, because traditional typing and amplicon-based sequencing methods may differentiate clonal complexes but miss gene content and single-nucleotide variants (SNVs) (14, 15). In contrast, shotgun metagenomics provides resolution of microbial communities at the strain and SNV levels (24). We used our previously validated strain tracking approach to identify strains of S. aureus and S. epidermidis present on our AD patients (16, 21). For S. aureus strain tracking, microbial reads were mapped against a database composed of 215 S. aureus genomes, of which 61 representatives are shown in Fig. 3A. Fig. 3. S. aureus–predominant individuals are often colonized with a single S. aureus strain. (A) Dendogram of 61 representative S. aureus strains based on SNVs in the core genome. Strains labeled in red were isolated from patients in (B). Colored blocks correspond to genomes of the same clade. Phylogenetically distant clade F1 is shown as an outgroup because it was recently reclassified as Staphylococcus argenteus (34). (B) For S. aureus–predominant individuals (N = 5), S. aureus clade relative abundances in bilateral Acs and Pcs for AD disease states, flare and post-flare. Colors correspond to those in (A). (C) For combined samples of all sites/time points of individuals in (B), bar charts show the number of SNVs per individual that are mono-, bi-, and triallelic. (D) Venn diagram showing the number of genes shared between isolates from patients in (B), indicated in red in (A). In contrast to the heterogenous communities of Propionibacterium acnes and S. epidermidis strains observed in healthy adult skin (16, 21), the more severe AD patients were markedly colonized with a single clade of S. aureus during disease flares (Fig. 3B, fig. S9, and table S8). In four of the five severe AD flare patients, this colonization with a single strain persisted in the post-flare but at notably lower mean relative abundances. AD patient AD11 was the exception, colonized by three different clades of S. aureus (E17, E7, and B1), with only clade E17 predominating during a flare. Notably, these more severe AD patients were colonized with distinct S. aureus clades. This supports previous studies demonstrating that AD patients do not share a single dominant S. aureus clone (11, 26–28). The variation in the clonal S. aureus clades colonizing AD patients raises the possibility that this heterogeneity may contribute to the differential course and/or therapeutic responses of AD patients. To confirm our strain tracking results, we used a complementary approach in which SNVs were identified in the S. aureus core genome (1.9 Mbp shared between all sequenced S. aureus). To power this analysis, we combined all sites and time points for each patient. In total, we identified 38,867 variant positions in the S. aureus core or ~10,000 single-nucleotide polymorphisms per patient. We then used the degree of polyallelism in each individual to infer genetic heterogeneity or the presence of multiple S. aureus strains. We calculated the number of mono-, bi-, and triallelic SNVs for each patient (Fig. 3C). Consistent with the strain tracking results, SNVs in clonal S. aureus–colonized AD patients were monoallelic at 93% of sites, whereas heterogeneous patient AD11’s SNVs were monoallelic at only 53% of sites. S. aureus isolates cultured from each of the more severe AD flare patients underwent whole-genome sequencing to confirm that the cultured patient isolates grouped into the respective clades predicted by the strain tracking of the metagenomic data (Fig. 3A). Colony picking from cultured swabs for patient AD11 isolated a representative from the dominant E17 and the nondominant B1 clades. On the basis of standard sensitivity testing methods and whole-genome sequencing analysis, five of the six S. aureus isolates from the more severe AD patients were methicillin-sensitive S. aureus (MSSA), consistent with higher incidences of MSSA than methicillin-resistant S. aureus (MRSA) cultivated from AD patient skin (29–31). Comparative genomic analysis of these six S. aureus strains revealed extensive heterogeneity in the gene content, as predicted on the basis of the initial mapping of the shotgun metagenomics sequences to disparate phylogenetic clades. The genome of a single S. aureus isolate encodes ~2500 genes, of which ~85% (2128 genes) are present in every strain’s genome and constitute the functional core (Fig. 3D); the remaining ~300 genes derive from the flexible pangenome composed of 1020 genes. We looked for functional enrichment in noncore versus core genes to identify pathways that were the most variable between our isolates (table S9). In doing so, we identified the KEGG pathways ko05150 S. aureus infection, ko00906 carotenoid biosynthesis, and ko01501 β-lactam resistance as functionally enriched in the variable component of the pangenome. With a targeted search, enterotoxin genes, previously shown to exacerbate AD (32), were differentially present in the six AD patient strains of S. aureus; variable presence of toxin genes is a trend consistent across S. aureus genomes (33). The five genes in the carotenoid biosynthesis pathway were present in all genomes but AD01.F1; this isolate is the most closely related to strain MSHR1132 that was recently reclassified as S. argenteus and can be visually distinguished by its white pigment versus yellow pigment (34). Finally, variability of genes in the β-lactam resistance family, including the mec cassette, was consistent with our previous result that only isolate AD11.E17 was an MRSA. Overall, this strain-level gene variation generates additional questions regarding the potential role of specific strains on disease pathogenesis and host factors on clonal strain selection.

Heterogeneous S. epidermidis strain communities are detected in AD and controls To further address the microbial community structure, we explored whether AD patients harbored heterogenous communities of S. epidermidis on skin. For S. epidermidis strain tracking, microbial reads were mapped against a database composed of 61 sequenced, phylogenetically diverse S. epidermidis genomes (Fig. 4A). As seen with healthy adults (21) and children, AD patients’ S. epidermidis communities at both flares and post-flares were composed of multiple different strains from diverse clades of the phylogenetic tree (Fig. 4B, fig. S10, and table S10). This directly contrasts with the identification of clonal S. aureus communities. This heterogenous S. epidermidis strain diversity was observed for both the more severe and less severe AD flare patients (fig. S10). However, analysis of the S. epidermidis strain composition in this cohort and our previous cohort of healthy adults (21) revealed a clustering of the less severe AD flare patients (Fig. 4C). Specifically, both unsupervised clustering and principal coordinate analyses identified S. epidermidis clades A29 and A30 as contributing to the clustering of the less severe AD patients and clade A20 as contributing to the clustering of the healthy adults (Fig. 4D). In contrast, the S. epidermidis strain diversity in healthy control children and more severe flare patients was intermixed. Fig. 4. S. epidermidis–predominant individuals are colonized by a heterogenous community of S. epidermidis strains. (A) Dendogram of S. epidermidis strains based on SNVs in the core genome. Strains isolated from patients in our study are labeled in red. Similar colors represent closely related strains that were grouped into 14 clades. Starred (*) isolates are nosocomial in origin. (B) For S. epidermidis–predominant individuals (n = 6), S. epidermidis strain relative abundances in AcPc for AD disease states, flare and post-flare. Colors correspond to those in (A). (C) Heatmap shows mean relative abundance of each clade across all sites in S. aureus– and S. epidermidis–predominant AD patients, healthy adults, and healthy children. (D) In principal component (PC) analysis, clades A20, A29, and A30 drive separation between S. epidermidis–predominant AD patients and healthy adults. S. epidermidis clades A29 and A30 were enriched in strains originally collected from nosocomial infections rather than as skin commensals (indicated with asterisks in Fig. 4A) (35). Comparative genomic analysis of nosocomial isolates and the other strains revealed higher relative abundances of the SCCmec cassette (35), which encodes genes necessary for methicillin resistance, in the nosocomial isolates. To further evaluate the S. epidermidis strains in this cohort, isolates were cultured from swabs collected from less severe patients AD05 and AD10. Whole-genome sequencing confirmed the patient isolates as members of the A29 and A30 clades, respectively (Fig. 4A, red). Consistent with the trend of increased antibiotic resistance genes observed through genomic analysis, these patient isolates were methicillin-resistant S. epidermidis. A potential explanation for the overrepresentation of isolates genomically similar to nosocomial strains in less severe AD flare patients may be that these S. epidermidis strains outcompete commensals and/or S. aureus in inflammatory or non–steady-state conditions or that antibiotic usage in these patients may have selected for antibiotic resistance genes.