Research ArticleMicrobiome

Staphylococcus aureus and Staphylococcus epidermidis strain diversity underlying pediatric atopic dermatitis

See allHide authors and affiliations

Science Translational Medicine  05 Jul 2017:
Vol. 9, Issue 397, eaal4651
DOI: 10.1126/scitranslmed.aal4651

Diving into the depths of atopic dermatitis

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 (TH2) and TH17 cells. Integrating high-resolution sequencing, culturing, and animal models demonstrated how functional differences of staphylococcal strains may contribute to the complexity of AD disease.


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).

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 (811), 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).

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.


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 (1719). 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 (68, 1113). 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, 2628). 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 (2931).

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.

Strains elicit differential cutaneous immune responses in a murine model

Whereas S. aureus has been tightly linked with AD, it is still debated whether S. aureus is a cause or effect, that is, whether S. aureus can elicit and/or worsen AD skin disease or is a bystander that flourishes with increased access to extracellular matrix or other products of inflammation in eczematous skin (36, 37). By observing that individual strains of S. aureus predominated during AD flares in our more severe flare patients, we sought to investigate whether these clonal strains elicited a biological response distinct from other strains of staphylococci. By harnessing the combined power of shotgun metagenomic sequencing of clinical samples and whole-genome sequencing of bacteria cultivated from concurrently collected skin swabs, we next analyzed (i) whether strains associated with AD flares would be sufficient to elicit skin inflammation in the absence of any known genetic predisposition or prior barrier disruption and (ii) whether there were strain-specific differences. To test this, we topically applied staphylococcal strains cultivated from AD patients and healthy controls onto intact skin of specific pathogen–free C57BL/6 wild-type mice, with a method previously developed to test the immune response to skin commensals. We individually tested 10 phylogenetically distinct S. aureus isolates: 6 S. aureus isolates cultivated directly from the flared skin of patients with more severe flares, 2 S. aureus isolates from the skin of the less severe patient AD07’s flare time point, 1 S. aureus isolate from a healthy control, and 1 common pathogenic S. aureus USA300 FPR3757 isolate (highlighted in red in fig. S9A). In addition, we tested three S. epidermidis isolates from AD patients: a representative each from clades A29 and A30, which predominated in the skin of less severe AD patients, and a representative from the ubiquitous B clade (highlighted in red in fig. S10A). In contrast to the noninflammatory responses observed after association with either skin commensals (38, 39) or AD patient S. epidermidis isolates, topical application of the S. aureus isolates, particularly those associated with more severe AD flare patients, was sufficient to induce epidermal thickening and inflammatory responses (Fig. 5, A and B, and fig. S11A) as well as immune cell infiltrate composed of neutrophils and eosinophils (Fig. 5C and fig. S11, B and C). The USA300 isolate, commonly used as a representative S. aureus in functional experiments, induced only a modest immune response as compared to many of the isolates cultivated from severe AD patients, underscoring the importance of using matched clinical isolates.

Fig. 5. Topical application of AD isolates induces AD-like cutaneous immune responses in murine models.

(A) Representative histological images of the ear pinnae of mice associated with tryptic soy broth (TSB); S. aureus AD04.E17, HC.B1, and USA300; or S. epidermidis A10.A30. Dotted line indicates separation between the epidermidis and dermis. Scale bars, 50 μm. (B) Epidermal thickness of ears after topical association of patient AD isolates. Color indicates origin and species of the isolate. (C) Absolute numbers of skin eosinophils, gated on Lineage, Ly6G, MHCII, CD64, and SiglecF+. (D) Absolute numbers of skin CD45+ TCRβ+ CD4+ cells. (E) Absolute numbers of skin IL-13+ TCRβ+ CD4+ cells. (F) Absolute numbers of skin IL-17A+ TCRβ+ CD4+cells. (G) Frequencies of IL-13+ and IL-17A+ CD4+ cells from mice in (B). Results are cumulative data from two or three independent experiments, with three mice per group. *P < 0.05, **P < 0.01, and ***P < 0.001, as calculated by analysis of variance (ANOVA) with multiple comparison correction.

In addition to innate immune cells, infiltration of T cell receptor (TCR) αβ+ and γδlow cells was also observed in mice colonized only with specific S. aureus strains (fig. S12A). Most of the TCRβ+ cells were CD4+ with variable effector potential, depending on the associated isolate (Fig. 5D). Notably, four S. aureus isolates from more severe AD flare patients induced production of the cytokine interleukin-13 (IL-13) (Fig. 5E), which is commonly associated with allergic inflammation. Cutaneous T helper 17 (TH17) cells were also identified when mice were colonized with these four IL-13–inducing strains, in addition to AD07.B2 and USA300 (Fig. 5, F and G). Recent reports have identified the presence of TH17 cells in AD lesions (40, 41), particularly in Asian patient populations (42).

Similar to CD4+ T cells, the γδ T cells of mice associated with specific strains of S. aureus isolates also had the potential to make higher levels of IL-17A (fig. S12B). Notably, four of the S. aureus isolates [two from the more severe flare patient AD11 (AD11.B1 and AD11.E17), one from the less severe flare patient AD07 (AD07.E7), and one from a healthy child (HC.B1)] induced minimal immune responses in all categories. Overall, association of S. aureus strains isolated from more severe AD flare patients to wild-type mice without prior barrier disruption induced immune responses in the skin that were significantly greater than those induced with S. epidermidis or S. aureus isolates from less severe AD flare patients or controls. Thus, these findings suggest that specific strains of S. aureus may be sufficient to elicit and/or exacerbate skin inflammation as part of AD disease pathogenesis.


AD is a complex disease with many contributing factors, including skin barrier integrity, innate and adaptive immunity, and the microbiome. The heterogeneity of the course, severity, and clinical response in AD patients underscores the diversity of phenotypic presentations, as well as the probable differences in disease pathogenesis, within this one diagnosis. In addition to the various genetic susceptibility loci for AD, deeper investigation into the skin microbiome could provide a better understanding of the microbial heterogeneity of AD and its potential contributions to disease.

Although there have been many efforts to identify bacteria in AD skin, studies have generally relied on methods that do not distinguish microbes beyond the species level or can misclassify genomically distinct clones (14, 15). Here, we combined shotgun metagenomic sequencing of clinical samples with whole-genome sequencing of patient-derived isolates to investigate the microbial communities of AD skin down to the strain- and SNV-level resolution. Because topical anti-inflammatory and antimicrobial treatments alter the skin microbiota (6, 7), the baseline and flare time points in this cohort were strictly defined by skin preparatory regimens to capture the natural history of the skin disease and to avoid potential confounders. As compared to healthy controls, the AD patients exhibited marked skin bacterial dysbiosis during flares. This dysbiosis was related to the increased relative abundance of staphylococci, consistent with prior cohorts. On the basis of the disease severity (defined by objective SCORAD) during flares, we observed a strong correlation between severe AD flares and S. aureus relative abundances. These findings demonstrated that despite the relatively small numbers of subjects in this study, our cohort of patients is representative of other published patient cohorts as defined by validated diagnostic criteria.

Shotgun metagenomic sequencing enables strain-level examination of microbes within the broad microbial community of bacteria, fungi, and viruses. Strain tracking identified marked outgrowth of clonal S. aureus strains in the skin of flaring AD patients with more severe disease; these same strains persisted post-flare at lower relative abundances. Other methods have examined whether S. aureus expansion in the skin of AD flares was related to either proportional increases in the entire community of S. aureus strains or the increase of a single or a few dominant clones; however, these studies were limited by the inability to examine these possibilities in the context of the whole skin microbial community. Although the fungal and viral communities were not significantly different in this study, expansion of reference databases/genomes and studies into the microbial “dark matter” in metagenomic data may provide further insights into AD microbiota. Our findings demonstrate that AD skin flares in patients with more severe disease are tightly linked with clonal S. aureus isolates.

In addition to characterizing strain communities during the course of AD, we found that less severe AD patients were colonized with more methicillin-resistant strains, whereas the more severe AD patients were primarily colonized with methicillin-sensitive strains. Although methicillin resistance is not as common in AD as would be predicted on the basis of the high rates of S. aureus colonization in this disease, the finding of MSSA and methicillin-resistant S. epidermidis predominance may contribute to differential responses to therapies in AD patients (43). The contrasts between S. aureus and S. epidermidis observed in this study likely also relate to the differences in microbial genetics and population dynamics at both the species and strain levels. Additional investigations of these microbiome phenotypic differences may improve the understanding of AD pathogenesis and lead to more targeted therapeutics, including the potential use of commensals to protect against S. aureus (44). Birth cohort studies may address whether these patients acquired bacterial strains from family members and/or environmental sources as part of microbial inheritance (45). Testing of S. aureus strains in gnotobiotic mice, similar to Bacteroides gut commensal studies, may functionally address whether colonization by clonal S. aureus occurs through limited exposure or colonization resistance (46).

Using strains isolated from the skin of AD flares and a healthy control as well as a known laboratory strain, we examined the potential biological differences between staphylococcal strains. In a murine model without prior skin barrier disruption and with intact immunity, S. aureus strains from flare time points in more severe AD patients were sufficient to induce manifestations of skin inflammation, such as epidermal thickening and cutaneous infiltration of TH2 and TH17 cells. The magnitude of different immunologic effects varied depending on the isolated strain but was not strictly related to the disease severity of the source patient. Notably, murine colonization with either isolate AD11.B1 or AD11.E17 induced minimal immune responses, although patient AD11 had an objective SCORAD of 51.4. However, AD11 is also heterozygous for a null mutation in the FLG gene (S757X), suggesting that AD11’s strains of S. aureus may be immunogenic in the setting of an impaired skin barrier, which, as shown by previous studies, allows S. aureus to breach the epidermis into the dermis where it can trigger expression of proinflammatory cytokines (47). Caveats of these findings in the murine model are the relatively small number of isolates from this cohort that were fully sequenced and studied in the murine model and the observation of varied host responses when testing isolates from the same clade (AD04 and AD11), highlighting the need to examine a larger number of isolates including strains from similar and different clades and from healthy individuals and AD patients. An important additional limitation is the recognition that this murine model and others do not recapitulate the multiple complexities of human AD.

In mouse models, S. aureus enterotoxins have been shown to act as superantigens that can initiate TH17 responses (48), whereas S. aureus δ-toxin can induce degranulation of mast cells (49). These genes were both present in the noninflammatory S. aureus isolates, indicating that strain variability exists not only in gene content but also in gene expression. Because healthy control–associated S. aureus strains were limited in our cohort because of the small percentage of healthy individuals colonized with S. aureus, future studies with additional S. aureus isolates from healthy individuals are necessary to tease apart the mechanisms underlying functional differences between S. aureus strains. In the context of prior studies demonstrating cutaneous immunologic responses to skin commensals (38, 39) and exacerbation of eczematous skin in AD mouse models by S. aureus (38, 39, 49, 50), our findings demonstrate that staphylococcal strains may play an important role in AD disease progression in a strain-specific manner.

Here, we used shotgun metagenomic sequencing to examine strain-level microbial compositions of AD skin, coupled with whole-genome sequencing of patient isolates. With increasing recognition of highly individualized skin microbiomes (16), the presence of patient-specific strains underscores the individuality of the disease course and therapeutic response and may represent an opportunity for precision medicine. Our functional studies with cutaneous colonization of AD patient–associated strains of S. aureus and S. epidermidis demonstrated strain-specific differences in the ability to elicit histologic and immunologic alterations. AD typically has an age of onset in the first year of life when the human immune system is developing and being tuned by the endogenous microbial community. Recent studies have shown that early exposures can modulate host immunity to subsequent exposure and induce tolerance (51, 52). Thus, in light of the known links between severe AD and subsequent development of asthma and hay fever (“the atopic march”), targeted modulation of an AD patient’s particular staphylococcal strains has the potential to ameliorate the broader development of atopic disorders.


Study design

AD patients and similarly aged healthy controls were recruited to participate in a natural history study approved by the institutional review board of the National Human Genome Research Institute (NHGRI) ( Written informed consent was obtained from parents or guardians of all participating children. Patients were diagnosed with AD on the basis of the UK Working Party definition (53). Eligibility criteria included ages of 2 to 18 years, moderate to severe disease (objective SCORAD, ≥15), presence of ≥1 affected Ac or Pc at enrollment, and >3 weeks off of systemic antibiotics and corticosteroids (1719). After skin preparation regimen, standardized skin sampling was performed from prespecified skin sites bilaterally and at defined time points (baseline, flare, and post-flare). Skin samples for metagenomic sequencing and negative controls were obtained as previously described (16, 21), with additional swabs of the Ac, retroauricular crease, and the nares collected concurrently for subsequent culture analyses.

Microbiome sequencing and analysis

Procedures for library generation with Nextera DNA Library Prep Kit and sequencing 2 × 125–bp reads with a target of 15 million to 50 million clusters on an Illumina HiSeq instrument were performed as described previously (21). In total, for 18 individuals (11 patients and 7 controls) sampled at seven body sites at different stages of disease for AD patients, we obtained 422 samples and 2.26 trillion reads (or 191 Gbp) of nonhuman, quality-filtered reads.

Microbial reads were assigned taxonomic classifications as previously described (21). Included in the microbial reference genome database are 2342 bacterial, 389 fungal, 1375 viral, and 67 archaeal genomes. In addition, a staphylococcus database was compiled from 315 complete and draft genomes from the National Center for Biological Information (NCBI) ( as of October 2014. Nonhuman reads were separately mapped to both genome collections using Bowtie2’s “--very-sensitive” parameter -k 10 to retrieve the top 10 hits (54). The resulting alignment files were processed with Pathoscope v1.0 (55) to assign multiply mapped reads to their most likely genome of origin. Read hit counts were then normalized by genome and scaled to sum to 1. Coverages of each output genome were calculated using genomeCoverageBed in the bedtools suite (56). To reduce the effects of spurious classifications from low-abundance organisms, only species with ≥1% coverage of the genome were considered (21).

Strain tracking of S. aureus and S. epidermidis was performed as previously described (21). Briefly, reference databases for S. aureus and S. epidermidis were compiled from all complete and draft genomes available on NCBI, 215 and 61, respectively. For both species, whole-genome alignment, with nucmer (57), was then used to identify the “core” region shared between all sequenced strains. SNVs identified in these core regions were subsequently used to generate dendograms with PhyML 3.0. For strain tracking to avoid noise from other staphylococcal species, metagenomic reads were first filtered against the staphylococcus database minus the species being strain-tracked (Bowtie2: --very-sensitive, -score-min L,-0.6,0.006). The remaining reads were then mapped to each species database with Bowtie2 (--very-sensitive, -score-min L,-0.6,0.006, -k number of genomes) (54) with zero tolerance for mismatches. The resulting alignment file was then processed with Pathoscope (-theta_prior 10 x 10^88) (55) to deconvolute multiple mapping reads.

For SNV analysis as described previously (16), metagenomic reads were mapped against the S. aureus core genome using Bowtie2 (--very-sensitive). The resulting alignment file was sorted by SAMtools and then processed with GATK’s IndelRealigner (58). The corrected alignment file was analyzed with SAMtools and BCFtools to identify possible variants (samtools mpileup -uD -q30 -Q30; bcftools view -Abvcg, varFilter -D99percentileofcoverage -d4 -1 .00001 -4 .00001). Custom scripts were then used to filter possible variants on the basis of criteria described in (59).

For staphylococcal isolates, Nextera libraries were generated from the genomic DNA and sequenced using a paired-end 300-base dual index run on an Illumina MiSeq instrument to generate 1 million to 2 million read pairs per library for ~80× genome coverage. Reads for each isolate were assembled with MaSuRCA version 2.2.1 (60), SPAdes version 3.6.0, or SPAdes version 3.6.0 (61) plus Pilon version 1.13 (62) correction. Best k-mer length estimates on paired-end reads were evaluated using KmerGenie (version 1.6300) (63) and used in running the MaSuRCA assembler for each genome. For comparative genomic analysis, genome annotation was done using the Institute for Genome Sciences (IGS) Analysis Engine ( (64).

To identify the functional capacity of S. aureus isolates, we followed a procedure similar to that in (16). The IGS Analysis Engine was used for structural and functional annotation of the sequences ( (64). Manatee was used to view annotations ( Protein sequences were then clustered into nonredundant orthologs with USEARCH (-cluster_fast -id 0.50 -centroids) (65). These gene clusters were then annotated by BLASTp against the KEGG database. Distribution of genes between the S. aureus assemblies was visualized with jvenn (66).

Murine topical association with bacterial isolates and flow cytometric analyses

Experiments were performed with 6- to 12-week-old female C57BL/6 specific pathogen–free mice under an animal study proposal approved by the NHGRI Animal Care and Use Committee. Topical association of mice was based on (38, 39). Ten S. aureus strains (AD04.E17, AD06.E13, AD03.A2, AD01.F1, AD11.B1, AD11.E17, AD07.B2, HC.B1, USA300, and FPR3757) and 3 S. epidermidis strains (AD05.A29, AD10.A30, and AD01.B), isolated from AD patients, controls (HC.B1), or the American Type Culture Collection (USA300), were cultured in tryptic soy broth at 37°C for 18 hours and normalized using optical density at 600 nm (OD600) to achieve similar bacterial density (about 108 colony-forming units/ml). For topical association, a sterile epicenter Catch-All swab was moistened in liquid culture of the bacteria and then rubbed against the ears of mice until they became visibly moist.

Cells from the ear pinnae of mice were isolated as previously described (38, 39). Murine single-cell suspensions were incubated with fluorochrome-conjugated antibodies against surface markers CD4 (clone RM4-5), CD8β (eBioH35-17.2), CD11b (M1/70), CD11c (N418 or HL3), CD19 (6D5), CD45.2 (104), CD49b (DX5), CD64 (X54-5/7.1), Ly6G (1A8), MHCII (M5/114.15.2), NK1.1 (PK136), TCRγδ (GL3), TCRβ (H57-597), and/or SiglecF (E50-2440). All antibodies were purchased from eBioscience, BioLegend, or BD Biosciences. For detection of basal cytokine potential, single-cell suspensions from ear tissue were directly cultured ex vivo in a 96-well U-bottom plate. Cell acquisition was performed on a Fortessa flow cytometer using FACSDiva software (BD Biosciences), and data were analyzed using FlowJo software (Tree Star).


All statistical analyses were performed in R, and most of the graphs were generated with ggplot2 (67). Data are means ± SEM unless otherwise indicated. For all box plots, center lines represent the median, and edges represent the first and third quartiles. The nonparametric Wilcoxon rank-sum test was used to determine statistically significant differences between populations (wilcox.test in R). Where indicated, within-subject analysis was performed with option “paired=T” in wilcox.test. All P values were adjusted using p.adjust in R using Bonferroni (no. of comparisons, ≤10) or false discovery rate (no. of comparisons, >10) corrections. Statistical significance was ascribed to an α level of the adjusted P values of ≤0.05. Similarity between samples was assessed using the Yue-Clayton θ, which assesses the similarity between two samples on the basis of (i) the number of features in common between two samples and (ii) their relative abundances, with θ = 0 indicating totally dissimilar communities and θ = 1 indicating identical communities (68). For functional experiments, statistical significance was determined by ANOVA with multiple comparison corrections (aov and TukeyHSD in R).


Materials and Methods

Fig. S1. Seven sites sampled bilaterally on pediatric AD patients and control children.

Fig. S2. Full multikingdom taxonomic classifications for AD patients and controls.

Fig. S3. Full Malassezia species classifications for AD patients and controls.

Fig. S4. Full eukaryotic virus classifications for AD patients and controls.

Fig. S5. Full bacterial taxonomic classifications for AD patients and controls.

Fig. S6. Relative abundance of staphylococcal species in relation to total bacterial population for all sites in AD patients and controls.

Fig. S7. Correlation of various staphylococcal species mean relative abundance and objective SCORAD for all sites of patients.

Fig. S8. Relative abundance of staphylococcal species for all sites in AD patients and controls.

Fig. S9. S. aureus clades for AD patients and controls.

Fig. S10. S. epidermidis clades for AD patients and controls.

Fig. S11. Histologic and cutaneous innate immune cell responses with AD isolate association in a murine model.

Fig. S12. CD45+ cutaneous immune responses with AD isolate association in a murine model.

Table S1. Subject Tanner stage and disease severity for samples used in this study.

Table S2. Clinical metadata for the subjects in this study.

Table S3. Metadata table for all samples in this study.

Table S4. Multikingdom relative abundances (coverage, >1% of the reference genome).

Table S5. Malassezia relative abundance.

Table S6. Bacteria relative abundances (coverage, >1% of the reference genome).

Table S7. Staphylococcus species relative abundances.

Table S8. S. aureus clade abundances.

Table S9. S. aureus pangenome analysis for isolates cultured and sequenced in this study.

Table S10. S. epidermidis strain relative abundances.

References (6975)


  1. Acknowledgments: We thank P. Thomas, M. Park, E. A. Kennedy, S. Phang, A. Pradhan, V. S. Pillai, I. Bozhenko [Harris Corporation employee working under contract with the National Cancer Institute (NCI)], and K. Beacht for the underlying efforts, and Segre lab, M. C. Udey, and K. Nagao for the helpful discussions. This study used the high-performance computational capabilities of the NIH Biowulf Linux cluster. IGS Analysis Engine at the University of Maryland School of Medicine provided structural and functional annotation of genomes. Funding: This study was supported by NHGRI, NCI, and National Institute of Allergy and Infectious Diseases Intramural Research Programs. Sequencing was funded by a grant from the NIH (4UH3AR057504-02). Author contributions: A.L.B., Y.B., J.A.S., and H.H.K. designed the study and drafted the manuscript. Sequencing was carried out by NISC. A.L.B. analyzed microbial sequence data. C.D., S.K.B.C., O.J.H., S.C., and W.-I.N. performed the experiments and analyses. Competing interests: The authors declare that they have no competing interests. Data and materials availability: The sequencing data for this study are linked to the NCBI BioProject 46333.
View Abstract

Navigate This Article