Study Design and Subjects

Supplementary Materials

Methods

Study Design and Subjects

This study was designed to investigate the skin microbiome in ADand was approved bythe Institutional Review Boards atNational Jewish Healthin Denver and UCLA in Los Angeles. All study participants underwent a detailed history survey, physical examination, and disease severity assessment.AD was diagnosed with standardcriteria based on the American Academy of Dermatology and The NIH/NIAID Atopic Dermatitis Research Network.1, 2Disease severity was assessed byusing the Rajka-Langeland scoring system,which rates extent, course, and itch intensity separately.3 Healthy subjects were defined as having no personal or family history ofatopic diseases and no personal history of chronic skin or systemic diseases.Exclusion criteria for all subjects included having a fever ≥ 38.5 ºC; having bathed or showered after midnightbefore the day of sampling; using creams/lotions at the sites where the specimens would be taken within 24 hours of sampling; receiving oral antibiotics,a bleach bath, topical prescription medications (including but not limited to Elidel, Protopic, topical corticosteroids, or topical antibiotics) at the sites where the swabs would be collected within seven days of skin sampling; systemic immunosuppressive drugs (including cyclosporine or oral steroids), total body phototherapy (e.g., ultraviolet light B, psoralen plus ultraviolet light A, tanning beds) within 20 days prior to sampling.

Allparticipants were enrolled at National Jewish Health. Written informed consent was obtained from all participants(parents or guardiansof participating children).The subjects were divided into three age groupsin the analyses: young children (age ≤ 12), teenagers (age 13-17), and adults (age ≥ 18). The divisions of age groups are in accord with the Tanner stages.4

Sample Collection, Genomic DNA Extraction, 16S rRNA Gene Amplification, and Sequencing

Two swab samples, one from lesional skin and one from non-lesionalskin, were collected from the volar surface of the forearm from each AD patient. A non-lesional site was defined as a clinically normal-appearing area that is approximately 5 cm away from the edge of the lesional site that is being swabbed. One skin swab samplewas collected from the volar forearmof each healthy subject. All sampleswerecollected froma 5cm X 5cm skin area. The swabswere placed into 2.0 ml locking tubes containing 300 μl ATL buffer (Qiagen, Inc., Valencia, CA), and were transported to the laboratory and stored at -80˚C.

Genomic DNA was extracted using QIAamp DNA micro kit (Qiagen, Inc., Valencia, CA) with modifications. Bead beating for optimal bacterial cell lysis was added.The hypervariable regions V1-V3of the 16S rRNA genes were amplifiedfrom purifiedgenomic DNA using primers 27F (5’-TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGAGAGTTTGATYMTGGCTCAG-3’) and 534R (5’-GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGCATTACCGCGGCTGCTGG-3’). PCR amplification was performed according to the protocol developed by the Human Microbiome Project.5Nextera indices were added to the purified amplicons based on the Illumina MiSeq specifications(Illumina, Inc., San Diego, CA).16S rRNA amplicon libraries were purified, quantified using qPCR and pooled for sequencing on Illumina MiSeq platform(Illumina, Inc., San Diego, CA). Paired-end reads of 251 bp were obtained for each sample.

Microbial Taxonomic Composition Analysis

A data cleaning process was applied to all sequence data prior to analysis.Briefly, low-quality bases with a Phred quality value lower than 20were trimmed off the read ends. The read pairs were removed if any of the two reads was trimmed to shorter than 200 bpor had ≥3% uncertain bases.The 16S rRNA sequences with paired reads were mapped against Greengenes(v13.5, non-redundant precalculated OTU references, 97_otus from PICRUSt).6, 7The alignments were performed using Bowtie2 to search sequences with the best hits to the references.8We required that the paired-end reads align to the same reference sequence in different strand directions and within 600 bp from each other.The annotation of 16S rRNA references wasrefined at the species level by cross-referencing the SILVA rDNA database.9The taxonomic composition of each sample was inferred at the phylum level, genus level, and species level. Some species in the same genus could not be distinguished due to their high similarityin the V1-V3 regionsof the 16S rRNA genes (≥ 99% nucleotide identity), and werecombined in the analysis,such as Streptococcus salivarius, Streptococcus thermophiles, and Streptococcus vestibularis. Taxawith ≥1% relative abundance in each sample were considered as present in the sample.

Analysis of the Microbial Community

The microbial community evenness and richnesswere measured by alpha diversity, and the similarity between individual microbial communities was measured by beta diversity. The evenness and richness (alpha diversity) of the microbial community were measured by Shannon index, calculated as follows:

wherepiis the proportion of theith operational taxonomic unit (OTU) in the sample, and Ris the total number of OTUs.

The similarity between individual microbial communities(betadiversity) was determined by weighted UniFracdistance,10 which measures the similarity between samples by incorporating both the relative abundance of the members of the microbial communities and phylogenetic distances between observed organisms. The betadiversity analysis, coupled with multivariatestatistical methodsincluding principalcoordinate analysis (PCoA) and Analysis of similarities(ANOSIM), can potentiallyidentify factorsexplaining differences between microbial communities.11PCoA is very similar to principal component analysis(PCA), with the main difference being that it uses metrics of community-wide similarity (distances between samples).ANOSIM is a widely used distribution-free method of multivariate data analysis,and is employed in this study to test whether the similarity within groups is statistically different from the similarity between groups.12 In the microbial community analysis of the skin samples, we separatedthe samples into two or more groupsbased on each categoryof the clinical metadata, and used ANOSIM to test whether there are significant differences between the groups in each category. Since ANOSIM is a non-parametric multivariate analysis, in this study the statistical significance was determined based on1,000 permutations.Analyses of the alpha diversity (Shannon index) and beta diversity (weighted UniFrac),Principal Coordinate Analysis (PCoA),and the sequencing depthassessment by the rarefaction analysis were calculated or performed usingQIIME.13Non-parametric multivariate analysis(ANOSIM)was conducted in Mothur14with default parameters.

We calculated the associations among the prevalent genera based on their relative abundances across all samples from the AD patients and the healthy controls. The genera were clustered based on the Pearson's correlations between genera.15 We identified three distinct clusters of microorganisms as shown in Figure 2b. To investigate the associations of skin microorganismswithage,we calculated the fold changeof the relative abundance of each microorganism between the two age groups in healthy controls. We found that the clusters of microorganismswere associated with specific age groups (young children versus adults-teenagers), with statistical significance testing on members of clusters.

Analysis of the Functional ProfilesofSkin Bacteria

To determine whether the skin bacteria identified in different age groups or disease state encode distinct functional profiles, the distribution of KEGG orthologous groups (KO genes)in the dominant skin bacteria was analyzed. KO genes are defined based on genome annotation as orthologous genes across multiple species.16Microorganisms that were present in any sample with ≥1% relative abundance and from the four representative genera (Staphylococcus, Streptococcus,PropionibacteriumandCorynebacterium) were included in the analysis. These organisms included AD-associatedS.aureus and S. epidermidis,childhood-associated Streptococcusspp. (Streptococcus pneumoniae,Streptococcus agalactiae, Streptococcus thermophilus, Streptococcus sanguinis,Streptococcus gordonii,Streptococcus dysgalactiaeand Streptococcus mitis),andadult-associated P.acnes and Corynebacteriumspp. (Corynebacterium jeikeium,Corynebacterium aurimucosum and Corynebacterium kroppenstedtii). When availablein HUMAnN,17multiple genomes of each species were included to cover the genetic diversity of the skin bacteria. In total, 46 genomes were analyzed;15 genomes from S. aureus, 2 from S. epidermidis, 24 from Streptococcusspp., 2 from P.acnes, and 3 from Corynebacteriumspp.(Supplementary Table E6). A KO gene wasconsidered present in an organism if it was identified in more than half of the genomes ofthe species. The bacterial species were clustered based on the presence/absence of KO genesin the genomes using hierarchical clustering in Cluster 3.0.18The distance measurement was based on the Pearson's correlation between genomes of species.

Statistical Analysis

Data are represented as mean ± s.e.m unless otherwise indicated. Student's t-test with two-tailed distribution was used in statistical testing of differences between samples from different groups. Paired Student's t-testwas used in comparisons between lesionaland non-lesional samples from the same individuals. The p-values were adjusted for multiple testing of microbial taxa with p.adjust in R using false discovery rate.19

References

1.Eichenfield LF, Hanifin JM, Luger TA, Stevens SR, Pride HB. Consensus conference on pediatric atopic dermatitis. J Am Acad Dermatol 2003; 49:1088-95.

2.Beck LA, Boguniewicz M, Hata T, Schneider LC, Hanifin J, Gallo R, et al. Phenotype of atopic dermatitis subjects with a history of eczema herpeticum. J Allergy Clin Immunol 2009; 124:260-9.

3.Rajka G, Langeland T. Grading of the severity of atopic dermatitis. Acta Derm Venereol Suppl (Stockh) 1989; 144:13-4.

4.Marshall WA, Tanner JM. Variations in pattern of pubertal changes in girls. Arch Dis Child 1969; 44:291-303.

5.The Human Microbiome Project Consortium. A framework for human microbiome research. Nature 2012; 486:215-21.

6.DeSantis TZ, Hugenholtz P, Larsen N, Rojas M, Brodie EL, Keller K, et al. Greengenes, a chimera-checked 16S rRNA gene database and workbench compatible with ARB. Appl Environ Microbiol 2006; 72:5069-72.

7.Langille MG, Zaneveld J, Caporaso JG, McDonald D, Knights D, Reyes JA, et al. Predictive functional profiling of microbial communities using 16S rRNA marker gene sequences. Nat Biotechnol 2013; 31:814-21.

8.Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nat Methods 2012; 9:357-9.

9.Quast C, Pruesse E, Yilmaz P, Gerken J, Schweer T, Yarza P, et al. The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Res 2013; 41:D590-6.

10.Lozupone C, Lladser ME, Knights D, Stombaugh J, Knight R. UniFrac: an effective distance metric for microbial community comparison. ISME J 2011; 5:169-72.

11.Lozupone C, Hamady M, Kelley ST, Knight R. Quantitative and qualitative beta diversity measures lead to different insights into factors that structure microbial communities. Appl Environ Microbiol 2007; 73:1576-85.

12.Clarke KR. Non-parametric multivariate analyses of changes in community structure. Aust. J. Ecol. 1993; 18:117-43.

13.Caporaso JG, Kuczynski J, Stombaugh J, Bittinger K, Bushman FD, Costello EK, et al. QIIME allows analysis of high-throughput community sequencing data. Nat Methods 2010; 7:335-6.

14.Schloss PD, Westcott SL, Ryabin T, Hall JR, Hartmann M, Hollister EB, et al. Introducing mothur: open-source, platform-independent, community-supported software for describing and comparing microbial communities. Appl Environ Microbiol 2009; 75:7537-41.

15.Eisen MB, Spellman PT, Brown PO, Botstein D. Cluster analysis and display of genome-wide expression patterns. Proc Natl Acad Sci U S A 1998; 95:14863-8.

16.Kanehisa M, Goto S, Kawashima S, Okuno Y, Hattori M. The KEGG resource for deciphering the genome. Nucleic Acids Res 2004; 32:D277-80.

17.Abubucker S, Segata N, Goll J, Schubert AM, Izard J, Cantarel BL, et al. Metabolic reconstruction for metagenomic data and its application to the human microbiome. PLoS Comput Biol 2012; 8:e1002358.

18.de Hoon MJ, Imoto S, Nolan J, Miyano S. Open source clustering software. Bioinformatics 2004; 20:1453-4.

19.Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. R. Stat. Soc. Series B 1995; 57:289–300.

ExtendedFigure Legends

Figure 1. The microbiome composition in healthy and AD patients.

a) The skin microbiome of the healthy controls is shown at the phylumlevel (upper panel) and at the genus level (lower panel). Seven bacterial phyla were present in at least two subjects and 20 genera were present in ≥20% of the samples of any age group. Samples were ordered by subject’s age at sampling, and were divided into three age groups as indicated in the color scale at the bottom.b)The representative skin microorganisms differ with age in the healthy skin. The relative abundances of skin species P. acnes, Corynebacterium spp.and S. epidermidis were significantly higher in adults than in young children. In contrast,Streptococcus spp. was significantly more abundant in young children than in adults. S.aureushad very low relative abundance in all age groups of healthy subjects. Data were visualized as box plots, which show the 25th, 50th and 75th percentiles of the data from bottom to top.The Tanner stagesand finerage divisionsreflecting human physical development are indicated.c)The microbiome composition of AD skin.The microbiome compositions of lesional (upper panel) and non-lesional (lower panel)skin are shown at the genus level. The same 20 prevalent genera identified in the healthy controls (lower panel in Figure 1a) were detected in AD patients. Samples collected from the lesional skin and non-lesional skin of the same AD patients are shown in the same order. The samples from young children are shown on the left, and those from adults-teenagers are shown on the right.

Figure 2. The microbiome is different between young children and adults in both health and AD with age-specific co-occurrence of skin microorganisms.

a)A principal coordinate analysis shows significant differences in the skin microbiome between young children and adults in healthy skin and in the non-lesional and lesional skin of AD patients. The skin microbiome in teenagers is in transition from young children to adults.b)Age-specific co-occurrence of skin microorganisms. Based on the abundance profiles,the skin microorganisms were clustered into adult-associated, childhood-associated, and AD-associated groups. The correlation between organismswasmeasured based on the Pearson's correlation and is indicated by a gradient from red(co-occurring) to green (co-exclusive), asshown in the heat map on the left.Fold changesof the genera in relative abundance between healthy young children and adults-teenagersare indicatedon the right.Positive numbers indicate that the generawere more abundant in young children, and negative numbers indicatethat the genera were more abundant in adults. Genera are labeled with asterisks if they were significantly different between the two age groups (p<0.05).

Supplementary Figure Legends

Supplementary Figure E1. Rarefaction curve of 16S rRNA gene sequence data. The rarefaction curve suggested the sufficient sequencing depths provided in this study. Data are representedas mean ±s.e.m.

Supplementary Figure E2.The microbiome composition in healthy skin. In the skin microbiome of the healthy controls, 15species were present in ≥20% of the samples of any age group. Samples were ordered by subject’s age.

Supplementary Figure E3.Co-occurrence ofStreptococcusandStaphylococcusspecies.At the species level, we were able to identify multiple species within the genera ofStreptococcusandStaphylococcus. Except for S.aureus,species from the same genus were positively correlated with each other in relative abundance. S.aureuswasnegatively correlatedwith other species,including those from the samegenus.

Supplementary Figure E4. Unique functional profilesofage-specific and AD-associated skin bacteria. The functionsencoded in 46genomes of the representativeskin microorganisms were annotated withKO genes.AD-associated (S.aureusand S. epidermidis),childhood-associated (Streptococcusspp.),and adult-associated (P.acnes and Corynebacteriumspp.) skin microorganisms have distinct functional profiles. The number of annotated genomes in each species is labeled in brackets after the species name.

Supplementary Tables

Table E1.Study participants in different age groups.

Table E2.The taxonomic compositions of the healthy skin microbiome at the genus level

Table E3. The taxonomic compositions of the healthy skin microbiome at the species level

Table E4. The taxonomic compositions of the AD skin microbiome at the genus level

Table E5.The taxonomic compositions of the AD skin microbiome at the species level

Table E6. The KEGG orthologous groups encoded in skin bacterial species

Table E7.Correlations between the skin microbiome and host phenotypes