Microbe-host interplay in atopic dermatitis and psoriasis

ARTICLE Microbe-host interplay in atopic dermatitis and psoriasis Nanna Fyhrquist et al. # Despite recent advances in understanding microbial diversity in skin homeostasis, the relevance of microbial dysbiosis in inflammatory disease is poorly understood. Here we perform a comparative analysis of skin microbial communities coupled to global patterns of cutaneous gene expression in patients with atopic dermatitis or psoriasis. The skin microbiota is analysed by 16S amplicon or whole genome sequencing and the skin transcriptome by microarrays, followed by integration of the data layers. We find that atopic dermatitis and psoriasis can be classified by distinct microbes, which differ from healthy volunteers microbiome composition. Atopic dermatitis is dominated by a single microbe (Staphylococcus aureus), and associated with a disease relevant host transcriptomic signature enriched for skin barrier function, tryptophan metabolism and immune activation. In contrast, psoriasis is characterized by co-occurring communities of microbes with weak associations with disease related gene expression. Our work provides a basis for biomarker discovery and targeted therapies in skin dysbiosis. https://doi.org/10.1038/s41467-019-12253-y OPEN Correspondence and requests for materials should be addressed to H.A. (email: [email protected]). # A full list of authors and their affiliations appears at the end of the paper. NATURE COMMUNICATIONS | (2019) 10:4703 | https://doi .org/10.1038/s41467-019-12253-y | www.nature.com/naturecommunications 1 1234567890():,; Interactions between commensal or pathogenic microbes and the hosts they colonize are central to the maintenance of homeostasis and the initiation of disease1. This rapidly advancing field is now starting to bear the fruits of interdisciplinary efforts but our understanding of microbe−host interactions is still limited2. Skin represents a primary tissue interface, where accessibility to both microbiome species and underlying tissue provides a unique opportunity for studies into host−microbiome interactions3in the initiation and maintenance of atopic/allergic or autoimmunetype inflammation4,5. Recent advances in analyzing microbial gene sequences in healthy skin has provided a comprehensive understanding of the classes of microbes and their diversity occupying distinct topographical niches6. Pioneering studies have expanded into shotgun sequencing approaches to further probe the taxonomic diversity and biogeography of microbes in small cohorts of healthy individuals7. Commensal skin microbes control adaptive skin immune homeostasis through interaction with specific subsets of antigen-presenting dendritic cells (DCs) and effector T-cell populations8,9, supporting their own survival, and protecting against the overgrowth of pathogens. Insights into the association of distinct microbial classes with inflammatory skin disease and their impact on the host genome are just emerging. In atopic dermatitis (AD), as a model of atopic/allergic inflammatory disease, commensal skin microbes are associated with disease flares10. AD-related dysbiosis is frequently characterized by the colonization by Staphylococcus aureus, and simultaneous loss of other, potentially beneficial species. S. aureus colonizes skin effectively, and expresses several virulence factors11 with proven roles in the pathogenesis of AD by studying effects in cell12 and animal models13,14.S. aureus is associated with severity of the disease, and may be the result of a combination of detrimental effects from S. aureus on the one hand, and the loss of beneficial effects from other members of the skin microbiota on the other hand. Therefore, the use of antimicrobial therapeutics that target S. aureus may not be the optimal choice as they may also wipe out beneficial species or strains, and break mutualistic interactions between the skin and its microbiota. Psoriasis (PSO) appears to induce physiological changes at the lesion site, selecting for a specific microbiota. The psoriasis (PSO)-associated skin microbiome displays disease-specific features that might have a diagnostic value15–17, but whether the differential microbiota has pathophysiologic significance remains undetermined. Here, we present a large-scale, comprehensive analysis of the microbiome and microbiome-associated host transcriptome in skin of healthy volunteers (HV), AD, and PSO patients, revealing two distinct patterns of host−microbe interactions in chronic skin inflammation. We report a significant increase in the abundance of S. aureus and loss of anaerobic species in AD, while PSO is characterized by the co-occurrence of multiple organisms, including Corynebacterium and Finegoldia species. In AD lesions, which are abundantly colonized by S. aureus, toxin production and metabolic reprogramming are highly over-represented microbial functions. The host, in turn, responds to the changes in the microbiota via altered expression of genes related to barrier function, metabolic reprogramming, antimicrobial defense mechanisms and T helper type 2 (T H 2) signaling. In PSO, our results suggest that members of Corynebacterium may play a regulatory role, which is reduced in disease. Results The skin microbiotas in AD and PSO are highly distinct. Skin swabs were collected using a standardized protocol (Supplementary Fig. 1a−c). A total of 3.36 million 16S rRNA gene reads were analyzed using QIIME v1.8.018 resulting in the identification of 17,725 operational taxonomic units (OTUs) at 99.3% identity level. After removal of rare OTUs, 3342 remained. Blasting sequences against the Greengenes 16S rRNA gene database revealed healthy skin microbiomes consistent with previous reports6. An overall analysis indicated clear differences between AD and HV (Fig. 1a, b, Supplementary Fig. 2), confirmed by nonmetric multidimensional scaling (NMDS) analysis of the 95 most abundant OTUs (Supplementary Fig. 3a, Supplementary Table 1). AD also showed a reduction in diversity (Supplementary Fig. 3b). A total of 51 highly abundant OTUs showed significant differences between HV, AD and PSO (Fig. 1a). We additionally carried out a careful analysis of the influence of confounding factors, including age, anatomical location, gender and clinical center (Supplementary Table 2, Supplementary Fig. 3c). After correction for the confounding effects, 8 out of 17 confounded OTUs from AD and 5 of 13 OTUs from PSO retained significance. For example, C. kroppenstedtii was associated with age, S. aureus with anatomical location and Lactobacillus sp. with gender, but they all remained significantly associated with disease after correction (p< 0.01). Associations were tested using the Kruskal−Wallis test for body site and institution, the Mann−Whitney Utest for gender, and Spearman correlation for age. The remaining 11 most significant OTUs are shown in Fig. 1c. The most significant result is an increase in the abundance of S. aureus in AD, associated with a significantly lower abundance of OTUs representing strictly anaerobic bacteria in AD (Fig. 1d). The loss of anaerobes in AD is not driven by S. aureus, indicated by repeated analysis of samples devoid of S. aureus (Supplementary Fig. 2d). Significant changes in PSO compared to HV include increases in the abundance of C. simulans and C. kroppenstedtii as well as Finegoldia and Neisseriaceae species. Lactobacilli, Burkholderia spp. and P. acnes were lower in abundance in both AD and PSO compared to healthy skin. Microbiota-based classification of AD and PSO. We next asked whether the skin microbiome discriminates inflammatory skin pathologies. A supervised learning pipeline was employed to construct classifiers to explore the key sets of microbial taxa. We identified 26 microbes discriminating AD vs. HV cohorts with an area under the curve (AUC) of 0.94 (class errors HV =0.03, AD =0.27, Fig. 2a). The most discriminative taxa were of the genus Staphylococcus, including S. aureus (Z=14.0), S. epidermidis (Z=5.8), Staphylococcus spp. (Z=−6.9) and Burkholderia spp. (Z=−7.5). An NMDS analysis highlights the role of S. aureus in differentiating AD (red dots) from HV (blue dots; Supplementary Fig. 4a). We identified 24 microbes that differentiate PSO vs. HV with an AUC of 0.85 (class errors HV =0.09, PSO =0.33, Fig. 2b). The top discriminating microbes were C. simulans (Z=15.5), Neisseriaceae g. spp. (Z=6.9), C. kroppenstedtii (Z=5.5), Lactobacillus spp. (Z=−8.4) and Lactobacillus iners (Z=−3.5). An NMDS analysis demonstrated a clear separation boundary between HV (green dots) and PSO (green dots; Supplementary Fig. 4b). We identified 15 microbes that differentiate between AD and PSO with AUC of 0.86 (class errors AD =0.29, PSO =0.13, Supplementary Fig. 4c). The top microbes were S. aureus (Z=14.0), S. epidermidis (Z=3.9), Finegoldia spp. (Z=−6.5), C.simulans (Z =−5.0), and C.kroppenstedtii (Z=−4.4). An NMDS analysis highlights the importance of S. aureus in AD (Supplementary Fig. 4d). Communities of microbes associate with AD or PSO.To understand the interactions between communities of microbes under different disease states, we used network principles to express co-occurrence relationships. We found distinct ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-12253-y 2NATURE COMMUNICATIONS | (2019) 10:4703 | https://doi.org/10.1038/s41467-019-12253-y | www.nature.com/naturecommunications differences between the community structures of microbes associated with AD and PSO. For microbes associated with AD, SparCC19 correlation between taxa resulted in 19 species with a correlation above a threshold of 0.2 ( p< 0.05, Fig. 2c). S. aureus negatively correlated with species including Corynebacterium spp., S.epidermidis,Tepidimonas spp. and Phyllobacterium spp. Of the 24 microbes identified as important for PSO classification (Fig. 2b), 16 species showed significant SparCC correlation (SparCC > 0.2, p< 0.05, Fig. 2d). The most discriminant taxa C. simulans and C. kroppenstedtii displayed positive correlations 851917 851668 610043 205025 4449324 74351 4473201 164003 4353642 4348347 4467218 4346894 819937 4354809 4480063 4369229 837884 4301457 1107940 4456068 820692 360483 4482598 4021335 378096 4303697 4306540 4350124 4349859 4309323 1081372 654307 1096610 755148 625320 14278 4474056 505749 1003210 912906 4327300 2110555 4349522 4349519 939571 370309 4476950 761594 279980 883806 496787 4294554 4471315 4421536 4408996 4318084 4299324 1004369 25259 511475 2901965 4446521 4481323 851925 4473664 370134 403853 4439089 25478 3208510 247720 1036883 4047452 4317476 4422405 103810 987144 4460228 4468125 4411187 1131523 285376 114999 940083 4327286 912997 362390 995817 4422718 13445 565753 3841245 282360 20360 4440643 Relative abundance 160 1 100% AD OTU # Healthy PSO Stat value 90% 80% 70% 60% 50% 40% 30% 20% 10% Other Proteobact._Alphaproteobacteria Proteobact._Betaproteobacteria Proteobact._Gammaproteobacteria Firmicutes_Clostridiales Firmicutes_Lactobacilales Firmicutes_Staphylococcus Bacteriodetes_Flavobacteriales Bacteriodetes_Bacteriodales Other Actinobacteria Actinobacteria_Propionibacteria Actinobacteria_Corynebacteria AD HV PSO AD HV PSO 02500 5000 AD HV PSO 0 200 400 AD PSO 0200 600 400 AD PSO 01020 AD PSO 0 50 100 150 AD PSO 0306090 AD PSO 020 40 60 Abundance [reads/sample] Abundance [reads/sample] Abundance [reads/sample] Abundance [reads/sample] Abundance [reads/sample] Abundance [reads/sample] Abundance [reads/sample] ** ** ** ** ** ** ** ** ** 7500 30 200 AD PSO 0200 400 Abundance [reads/sample] ** AD PSO 0 200 400 600 Abundance [reads/sample] ** AD PSO 0 50 100 150 Abundance [reads/sample] ** ** AD PSO 010 20 30 Abundance [reads/sample] ** 40 ** HV HV HV HV HV HV HV HV HV S. aureus Finegoldia sp. Lactobacillus sp. C. kroppenstedtii Anaerobe avg. relative abundance 566 786 766 86 6 a b c d Streptococcus sp. Corynebacterium sp. Corynebacterium simulans Corynebacterium kroppenstedtii Corynebacterium sp. Corynebacterium sp. Dermabacter sp. Micrococcus sp. Rothia dentocariosa Kocuria palustris Propionibacterium acne s Actinomyces Prevotella sp. Prevotella sp. Bradyrhizobium sp. Caulobacteraceae Paracoccus sp. Phyllobacterium sp. Enhydrobacter sp. Acinetobacter sp. Neisseriaceae Methylophilaceae Burkholderia sp. Variovorax paradoxus Ralstonia sp. Pelomonas sp. Finegoldia sp. Peptoniphilus sp. Anaerococcus sp. Anaerococcus sp. Anaerococcus sp. Anaerococcus sp. Anaerococcus sp. Blautia sp. Peptostreptococcus anaerobius Staphylococcus epidermidis Staphylococcus sp. Staphylococcus aureus Lactobacillus sp. Lactobacillus iners AD HV Anaerobe abundance in AD and HV =Actinobacteria =Proteobacteria =Firmicutes =Bacteroidetes =Cyanobacteria Burkholderia sp.C. simulans Neisseriaceae Staphylococcus sp. P. acnes Pelomonas sp. Anaerococcus sp. NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-12253-y ARTICLE NATURE COMMUNICATIONS | (2019) 10:4703 | https://doi .org/10.1038/s41467-019-12253-y | www.nature.com/naturecommunications 3 Fig. 1 Characterization of the skin microbiome in AD and PSO. The results show a typical range of skin microbiomes in HV (n=115) and significant changes in AD (n=82) and PSO (n=119). aAn evolutionary tree based on 16S rRNA gene sequences, abundance and statistical significance of the 95 most abundant OTUs. The blue color intensity in the heat map shows the relative abundance of each OTU. The three Staphylococcus OTUs (indicated by asterisks) were calculated using a wider scale, due to their relatively high abundance in certain samples. The length of the green vertical bars indicate nonparametric statistical score (Kruskal−Wallis test, FDR, p< 0.05). The color bars by each OTU number indicates bacterial phylum (green: Firmicutes; red: Proteobacteria; cyan: Bacteroidetes; blue: Cyanobacteria; orange: Actinobacteria). bThe most abundant bacterial groups depicted for HV, AD and PSO. cStatistical analysis (Mann−Whitney Utest (FDR, p< 0.05)) of the 11 OTUs showing the most significant changes in AD and/or PSO vs. HV after correction for confounding factors. The values on the xaxis are in number of reads, out of a total of 8495 reads/sample. The asterisks indicate statistically significant differences and correspond to p< 0.01 (**) and p< 0.001 (***). The center line in the boxplots corresponds to the median, the bounding box is the interquantile range (IQR) and the whiskers are defined as 1.5 times IQR. dStatistical analysis (Mann−Whitney Utest, p< 0.05) of the relative abundance of OTUs representing strictly anaerobic bacteria in AD lesions and HV. The x-axis units and the elements of the boxplots are as in (c). The source data files used to generate the present figure are available from the NCBI Sequence Read Archive under accession PRJNA554499 8 a b c d Positive correlation Negative correlation Log2 FC (disease-healthy) –4.5 Staphylococcus aureus Staphylococcus sp. Propionibacterium acnes Staphyloccous epidermidi s Corynebacterium sp. Corynebacterium sp. Bradyrhizobium sp. Pelomonas sp. Lactobacillus sp. Lactobacillus iners Comamonadaceae g. sp. Anaerococcus sp. Rothia dentocariosa Finegoldia sp. Caulobacteraceae g. sp. Tepidimonas sp. Phyllobacterium sp. Burkholderia sp. Caulobacteraceae g. sp. Actinomyces sp. Corynebacterium simulans Peptostreptococcus anaerobius Corynebacterium sp. Rothia mucilaginosa Burkholderia sp. Streptococcus sp. Rothia dentocariosa Neisseriaceae g. sp. Corynebacterium kroppenstedtii Corynebacterium kroppenstedtii Neisseriaceae g. sp. Lactobacillus sp. Lactobacillus sp. Methylophilaceae g. sp. Anaerococcus sp. Staphylococcus aureus Staphylococcus epidermidis Corynebacterium sp. Rothia dentocarious Fusobacterium sp. Corynebacterium durum Lactobacillus sp. Corynebacterium sp. Caulobacteraceae g. sp. Teipidimonas sp. Anaerococcus sp. Finegoldia sp. Phyllobacterium sp. Lactobacillus iners Corynbacterium sp. Paracoccus marcusii Propionibacteirum acnes Anaerococcus sp. Corynebacterium sp. Lactobacillues sp. Comamonadaceae g. sp. Pelomonas sp. Caulobacteraceae g. sp. Bradyrhizobium sp. Staphylococcus sp. Burkholderia sp. Corynebacterium simulans Neisseriaceae g. sp. Lepotrichia sp. Corynebacterium kroppenstedtii Corynebacterium kroppenstedtii Peptostreptococcus sp. Streptococcus sp. Rothia mucilaginosa Corynebacterium sp. Peptostreptococcus anaerobius Anaerococcus sp. Actinomyces sp Neisseriaceae g. sp. Coprococcus sp. Rothia dentocariosa Lachonspiraceae g. sp. Fusobacterium sp. Corynebacterium sp. Lactobacillus sp. Methylophilaceae g. sp. Burkholderia sp. Pelomonas sp. Lactobacillus iners Lactobacillus sp. Variable importance (Z) –10 Selection frequency, % 100 80 60 40 20 0 100 Variable importance (Z) 100 –10 Fig. 2 Microbiota-based classification of AD and PSO. aVariable importance for the best set of discriminatory AD (n=82) taxa identified through Random Forest feature selection and classification (RF) analysis. Bars are colored by selection frequency (Red =OTU selected in all folds, blue =OTU selected in one fold). bVariable importance for the best set of discriminatory PSO (n=119) taxa identified by RF analysis. cCo-occurrence network of AD-associated OTUs selected by RF selection. Pairwise correlations were calculated over AD samples using SparCC and OTUs with a mean Z> 0.2 are shown. A connection represents a correlation >0.2 and significant p< 0.05 correlation. Solid and dashed lines respectively represent positive and negative correlations. The size and color of each node is proportional to the log2 fold change between healthy and disease. dCorrelation network of PSO-associated OTUs selected by RF. The source data files used to generate the present figure are available from the NCBI Sequence Read Archive under accession PRJNA554499 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-12253-y 4NATURE COMMUNICATIONS | (2019) 10:4703 | https://doi.org/10.1038/s41467-019-12253-y | www.nature.com/naturecommunications with Streptococcus spp., P.anaerobius and, Anaerococcus spp., Neisseriaceae g. spp., and Rothia dentocariosa respectively. Comparison across microbial interactions in PSO indicate that rather than a single species dominating the microbial landscape (as in AD), we observe multiple species associated with this disease type. Overlaps and distinctions of AD and PSO transcriptomes. The transcriptomes of full-thickness skin biopsies were defined based on a stringent 10−5FDR level and a fold change (FCH) of 1.5, or greater. A principal component analysis (PCA) revealed a clear separation between HV, AD, and PSO (Fig. 3a). Differential gene expression analysis identified 1232 genes differing between AD Term z-score Term z-score Term z-score c PSO/AD PSO/HV AD/HV abAD/HV PSO/HV PSO/AD 120 (4.7%) 1206 (43.8%) 854 (31%) 197 (7.2%) 51 (1.9%) 268 (9.7%) 48 (1.7%) d HV AD PSO 0 5 10 15 IL36A p < 0.0001 p < 0.0001 HV AD PSO 2 4 6 8 10 KYNU p < 0.0001 p < 0.0001 HV AD PSO 0 5 10 15 PI3 p < 0.0001p < 0.0001 HV AD PSO 2 4 6 8 10 NOS2 p < 0.0001 HV AD PSO 2.5 3.0 3.5 4.0 4.5 5.0 IL17A p < 0.0001 HV AD PSO 2 4 6 8 CCL1 p < 0.0001 p < 0.0001 HV AD PSO 2 4 6 8 10 CCL18 p < 0.0001p < 0.0001 HV AD PSO 2 3 4 5 6 7 IL22 p < 0.0001 p < 0.01 HV AD PSO 5 6 7 8 9 10 IL34 p < 0.0001 p < 0.0001 AD_LES CTRL PSO_LES –log P-value –log P-value –log P-value Ratio Ratio Ratio 13 2 0.25 0.10 9 2 0.54 0.16 0.194 0.063 4.9 2 –2.4 3.9 –2.09 3.8 –2 3.3 Th2 pathway Th1 pathway Dendritic cell maturation iCOS-iCOSL signaling in T helper cells Role of PRR in recognition of bacteria and viruses CD28 signaling in T helper cells p38 MAPK signaling Phospholipase C signaling Interferon signaling TREM1 signaling PKCθ signaling in T lymphocytes IL-8 signaling MIF-mediated glucocorticoid regulation Role of NFAT in regulation of the immune response MIF regulation of innate immunity LXR/RXR activation iNOS aignaling Tec kinase signaling Death receptor signaling Inflammasome pathway Oncostatin M signaling Cell cycle: G2/M DNA damage checkpoint regulation Interferon signaling LPS/IL-1 mediated inhibition of RXR function LXR/RXR activation Inflammasome pathway Dendritic cell maturation Role of IL-17A in psoriasis PPAR signaling Eicosanoid signaling iNOS signaling Mitotic roles of polo-like kinase Estrogen-mediated S-phase entry Toll-like receptor signaling p38 MAPK signaling Phospholipase C signaling Glioblastoma multiforme signaling TREM1 signaling Role of NFAT in regulation of the immune response Interferon signaling LXR/RXR activation p38 MAPK signaling Toll-like receptor signaling Role of PRR in recognition of bacteria and viruses Cholecystokinin/gastrin-mediated signaling IL-6 signaling RU Log2 RU Log2 RU Log2 RU Log2 RU Log2 RU Log2 RU Log2 RU Log2 RU Log2 Second component — inertia explained: 9.8% First component — inertia explained: 38.4% NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-12253-y ARTICLE NATURE COMMUNICATIONS | (2019) 10:4703 | https://doi.org/10.1038/s41467-019-12253-y | www.nature.com/naturecommunications 5 and HV, and 2525 between PSO and HV, and 1051 genes shared between AD and PSO (Fig. 3b). For functional insight into the disease-specific genes, we performed Ingenuity Pathway Analysis (IPA), revealing significant overrepresentation of T H 2 and T H 1 signaling, dendritic cell maturation and iCOS-iCOSL signaling in helper T cells in AD, and Interferon signaling, LPS-IL-1-mediated inhibition of RXR function, the Inflammasome pathway and Th17 signaling in PSO. The levels of Interferon signaling and p38 MAPK signaling distinguished PSO from AD (Fig. 3c). TNF and IFNG were identified as upstream regulators in both AD and PSO, whereas IL4 and IL13 emerged as unique features in AD, underpinning the central role of T H 2 -associated signaling in AD (Supplementary Fig. 5b). Gene Ontology (GO) enrichment analysis highlighted chemotaxis, inflammatory response, and extracellular matrix organization in AD, and leukocyte activation in both AD and PSO (Supplementary Fig. 6). Among top upregulated genes in AD and PSO we observed inflammatory mediators (S100 proteins, defensins, matrix metalloproteinases, IL-1 family cytokines), T helper-related genes (CCL1,CCL18, IL17A, IL22, PI3/Elafin), barrier genes (KRT16,SERPINB4,KLK9, FLG2,LCE5A,CLDN8), and genes involved in tryptophan (trp) metabolism (KYNU). Top downregulated genes included IL34, anti-inflammatory IL37, and NOS2 (Fig. 3d, Supplementary Fig. 7). The S. aureus-induced host gene signature in AD. Taking advantage of our large microbiome and transcriptome datasets, we addressed the interplay between host and microbes. We stratified AD samples into “high”and “low”groups, based on microbial abundance, including top (n=27) and bottom (n=25) tertiles in the analysis. The low abundance samples were devoid of S. aureus, whereas high abundance samples exhibited in all cases a high abundance of S. aureus (87–99%). First, we explored the microbial gene repertoire of S. aureus high and low groups. Using PiCRUST that predicts functional content based on the 16S rRNA gene marker, we observed the enrichment of bacterial toxins, the two-component system and glycolysis in the S. aureus high group compared with the low group (Supplementary Fig. 8a−b). The majority of the functions were contributed by S. aureus OTUs (Supplementary Fig. 8c). To validate OTU-based predictions, we performed whole genome sequencing (WGS) on a limited sample set, showing close agreement between the two independent sequencing methods (Supplementary Fig. 9). WGS data confirmed the enrichment of major bacterial toxins, hld and plc (alphaand delta-toxin, respectively), and genes associated with galactose metabolism, the phosphotransferase system, the two-component system, and glycolysis and gluconeogenesis, in the S. aureus high samples (Supplementary Fig. 8c). Next, we investigated gene expression profiles in the underlying skin. Comparison of the transcriptomes between S. aureus high and low samples revealed a set of 256 significant genes (FDR < 0.05, FCH ≥1.5) (Supplementary Fig. 10a, Supplementary Table 3). To explore whether the S. aureus-regulated genes were relevant to global features of AD pathophysiology, we created a co-expression network based on pairwise Pearson correlation (r> 0.7), using all significant AD-associated genes. Network community detection20 identified ten distinct modules, enriched for disease-relevant functions (Fig.4a, Supplementary Fig. 10c). Projecting S. aureus-regulated genes onto the AD network revealed significant enrichment in genes that mapped to modules M1 and M5 (hypergeometric test, FDR < 0.05), associated with keratinocyte differentiation, and extracellular matrix organization, respectively (Fig. 4b). Functional analysis of the S. aureus-regulated genes using GO and IPA revealed the enrichment of keratinization and skin development (Fig. 4c), and T H 17 signaling and tryptophan (trp) degradation (Supplementary Fig. 11a), respectively. IPA predicted IL1B,TNF and IFNG as top upstream regulators (Fig. 4d) and leukocyte migration and development of epithelial tissue as downstream effects (Fig. 4e). While upregulated genes were enriched for inflammatory signaling (Supplementary Fig. 11b), downregulated genes were over-represented for mainly skin development (Supplementary Fig. 11c). Top genes included skin barrier and antimicrobial factors (S100A7,DEFB4A/B, S100A9, MMP12, FLG2, CLDN8, ADAM12), components of trp metabolism (KYNU, TDO2, KMO), immune activation (IL1B, CCL2, CCL19), and T H 2 signaling (IL4R, IL5, IL13, PI3,TNFRSF4, CCR4). Moreover, HIF1A and its targets HK2 and PFKP were among the significantly regulated genes (Fig. 4e, f, Supplementary Fig. 10b, d−e). Since the kynurenine pathway was enriched in the S. aureus “high”samples, we performed transcriptomic reconstruction of trp breakdown, indicating the accumulation of 3hydroxyanthralic acid (3-HAA) in the skin (Supplementary Fig. 12a). Further, we investigated to what extent AD-associated strains of S. aureus may depend on trp, and isolated 32S. aureus strains from AD patients with moderate to severe manifestation of the disease. We found that 66% of the isolated strains grew independent of trp (Supplementary Fig. 12b, Supplementary Table 4), and this observation was further supported by our WGS data (73% trp biosynthesis-related genes in the S. aureus high samples) (Supplementary Fig. 12c). Finally, to explore the impact of S. aureus on gene targets identified by this study, organotypic human epidermal equivalents were topically exposed to 106CFU S. aureus for 24 h, followed by measurement of gene expression by qPCR (Supplementary Fig. 13a−c). S. aureus significantly induced the expression of DEFB4,PI3,IL4R and S100A9 (Supplementary Fig. 13d). Integrative analysis of the PSO microbiome and transcriptome. Next, we combined the PSO microbiome with associated skin transcriptomes. Testing the top discriminating microbes in PSO in terms of their ability to partition the skin transcriptome did not yield significant results. Therefore, we constructed a coexpression network based on PSO differentially expressed genes, and partitioned it into 12 modules as described above20 (Supplementary Fig. 14a). We identified associations between the top 25 most differentially abundant taxa between PSO and HV, and module eigengenes derived from the PSO co-expression Fig. 3 The AD and PSO skin transcriptomes. aProjection of AD (red, n=82), PSO (orange, n=119) and HV (green, n=115) transcriptome profiles in the subspace spanned by the two first components of the principal component analysis (PCA) performed on the 1000 most variant genes. bVenn diagram of the differentially expressed genes in the AD vs. HV, PSO vs. HV and PSO vs. AD contrasts (identified using the Limma linear model and empirical Bayes method, cut-off: log2 FC > 0.58 and FDR p<1×10 −5, corrected using the Benjamini−Hochberg method). cIPA canonic pathway analysis of significantly enriched functions in AD and PSO gene signatures and the PSO vs. AD contrast.dStatistical analysis (unpaired ttest) of selected top upand downregulated genes. The center line in the dot plots corresponds to the mean, and the error bars to the standard deviation. The source data files used to generate the present figure are available from EBI ArrayExpress under accession E-MTAB-8149 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-12253-y 6NATURE COMMUNICATIONS | (2019) 10:4703 | https://doi.org/10.1038/s41467-019-12253-y | www.nature.com/naturecommunications network, using linear models implemented in the MaAsLin package21. A linear model controlling for body site, institution, age and gender effects was fit for each microbe−module eigengene pair. Overall, six associations were identified (FDR < 0.20) (Supplementary Fig. 14b). Negative associations were identified between Corynebacterium spp. and three co-expression modules; M2, M3 and M10. M3 was strongly associated with the cell cycle, whereas M2 and M10 were enriched for inflammatory pathways including interferon signaling, and T H 1 and T H 2 activation (Supplementary Fig. 14c). These results suggest that Corynebacterium spp. may play a regulatory role which is reduced in disease however. Overall, our results indicate that microbe−host associations in PSO are considerably less well defined than observed in AD. Microbiome associates with clinical severity in AD. To investigate the effect of microbiota and gene expression on clinical S. aureus low S. aureus high 8.5 9.0 9.5 10.0 10.5 11.0 HIF1A p < 0.0001 S. aureus low S. aureus high 3 4 5 6 7 ADAM12 p = 0.002 S. aureus low S. aureus high 4 6 8 10 12 14 16 S100A9 p = 0.0025 S. aureus low S. aureus high 3 4 5 6 7 8 9 KYNU Atopic dermatits FABP5 PER1 SERPINB3 CCBE1 JAK3 CCL19 SLIT3 ADRB2 CDH3 SERPINB1 PCOLCE2 SEMA3E IL21R ZC3H12A TIMP3 IL37 SELE NR3C2 REN CCL2 BMP2 PLAT GJB2 XDH HRH2 S100A8 S100A9 DEFB4A/DEFB4B IL4R LEPR Development of epithelial tissue MGP INHBB ANGPTL1 ANGPTL4 TNMD PGF MMP12 S100A7A p = 0.0001 S. aureus low S. aureus high 7 8 9 10 11 12 13 FLG2 p = 0.0054 S. aureus low S. aureus hi g h 6 7 8 9 10 11 IL4R p = 0.009 S. aureus low S. aureus hi g h 2.5 3.0 3.5 4.0 4.5 5.0 IL13 p = 0.0164 S. aureus low S. aureus hi g h 1.5 2.0 2.5 3.0 3.5 IL5 p = 0.002 S. aureus low S. aureus high 0.0 2.5 5.0 7.5 10.0 12.5 15.0 DEFB4A/DEFB4B p = 0.0037 GO biological process complete fEnrich –log p-value Keratinization (GO:0031424) Keratinocyte differentiation (GO:0030216) Epidermal cell differentiation (GO:0009913) Skin development (GO:0043588) Epidermis development (GO:0008544) Epithelium development (GO:0060429) System development (GO:0048731) Epithelial cell differentiation (GO:0030855) Tissue development (GO:0009888) Developmental process (GO:0032502) Animal organ development (GO:0048513) Anatomical structure development (GO:0048856) Multicellular organism development (GO:0007275) Multicellular organismal process (GO:0032501) Cellular developmental process (GO:0048869) Cell differentiation (GO:00301 Peptide cross-linking (GO:0018149) Biological_process (GO:0008150) Negative regulation of leukocyte apoptotic proc. (GO:2000107) Response to external stimulus (GO:0009605) 14.8 20.8 Upstream regulator Predicted activation state Activation z-score –log pvalue of overlap OSM Activated Cg IL1B Activated IL1A Activated NFKB1A TNF Activated EHF IFNG INHBA Activated IL17R CHUK DUSP1 AGER FOS ADAMTS12 Interferon alpha Activated IL1 IL17C HIC1 IL17A Activated S. aureus signature enrichment Top enriched GO term M1 1.07E-12 Keratinocyte differentiation M2 1 Immune response M3 1 Lipid metabolic process M4 1 Regulation of cell cycle M5 0.002 Extracellular matrix organization M6 0.645 Metal ion transport M7 1 Response to virus M8 1 n.s. M9 1 Regulation of type 2 immune resp. M10 1 Lipid metabolic process = S. aureus signature –1.6 12.93 a b cd e f RU Log2 RU Log2 RU Log2 RU Log2 RU Log2 RU Log2 RU Log2 RU Log2 RU Log2 4.5 NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-12253-y ARTICLE NATURE COMMUNICATIONS | (2019) 10:4703 | https://doi.org/10.1038/s41467-019-12253-y | www.nature.com/naturecommunications 7 severity, we performed feature selection coupled with multivariate regression. Genes and microbes were ranked according to their correlation with SCORAD or PASI indices and then regression models were trained onto the top Nfeatures (N=5–50 in increments of five), and the best performing model was selected. Prediction of SCORAD from microbial abundance was optimal with the top 35 species with Random Forest regression (Supplementary Fig. 15a). The accuracy was relatively modest (MAE = 12.35, correlation true vs. predicted =0.55) suggesting that variability in microbial abundance and clinical severity are likely explained by other factors. Amongst the most highly correlated species were Tepidimonas spp., Propionibacterium acnes, and S. aureus (Supplementary Fig. 15b). Most of the selected species were negatively correlated with the exception of S. aureus, indicating that microbial diversity may be inversely associated with clinical severity in AD. Diversity correlation with SCORAD revealed a weak but significant correlation (cor =−0.27, p= 0.01). Transcriptomic data outperformed microbial abundance for predicting SCORAD. The top model identified 15 genes using linear regression as the best set of predictive genes (Supplementary Fig. 15c), (MAE =9.84, correlation true vs. predicted = 0.66). Amongst the top genes were SEMA3D,IGSF10, LGR5 and CASP10 (Supplementary Fig. 15d). No associations between clinical severity and microbial abundance in PSO were identified. Finally, oozing and crusting in the AD lesions indicative of infection, correlated significantly with the abundance of S. aureus, supporting the link between colonization by S. aureus and clinical severity in AD (Supplementary Fig. 16). Discussion Microbe−host interactions may be central to skin homeostasis, and dysbiosis may drive disease10,15. However, to what extent the skin microbiota may associate with the host skin phenotype and underlying mechanisms, remains elusive. To achieve a better understanding of the dialogue between the skin and the skin microbiome, we used skin as a model and AD and PSO as proxies for T H 2-associated atopic and T H 17-associated autoimmune inflammation, respectively, giving us the opportunity to compare two types of inflammation and the associated microbiomes and transcriptomes. Here we present the findings of a large cohort (n=316) combining the analysis of skin microbiomes and associated transcriptomes, identifying unique gene profiles that characterize healthy vs. inflamed skin. We show that while AD is dominated by one single microbial species, multiple species associate with PSO, and the abundance of the dominating species in AD, S. aureus, correlates with disease relevant gene expression (overview of the main findings in Fig. 5). Analysis of the two patient groups shows distinct differences between microbiomes. The most significant change, the increase in S. aureus abundance in AD, is not present in all samples. A number of AD lesions contain little or no S. aureus, representing potentially different endotypes of the disease which deserve further investigation. The changes associated with psoriasis are more complex than for AD and involve many different bacteria, including Corynebacteria and Finegoldia. Through exploring classification models as a way to discriminate disease state, several of the OTUs were also strong predictors of disease status. Lactobacillus appeared to be consistently depleted in both diseases, whereas disease-specific species such as S. aureus in AD and C. simulans and C. kroppenstedtii in PSO were important for prediction in disease. Moreover, we identified two distinct community patterns by means of co-occurrence analysis. In AD, it was clear that S. aureus dominated the microbial landscape and negatively correlated with several skin commensals, such as S. epidermidis and Corynebacterium spp, and thus may be associated with the depletion of potentially regulatory or protective microbes. Importantly, recent studies have shown that S. epidermidis may specifically limit the growth of S. aureus22, and disease severity is inversely correlated with the abundance of S. epidermidis relative to S. aureus14. In contrast to the dynamics revealed in AD, the PSO-associated interaction network suggested that multiple, co-occurring species, is a more representative model. Previous studies have shown that while single taxa are unable to discriminate between PSO and healthy skin, the combined relative abundances of selected genera (Corynebacterium, Streptococcus and Staphylococcus) attain significance across the groups15. In our hands, the same genera display clear patterns of co-occurrence, increasing in relative abundance in PSO lesions compared with healthy skin. We observed a remarkable loss of strictly anaerobic bacteria in AD, indicating a switch in the skin microbiome from anaerobic to aerobic metabolism. Healthy skin is normally O 2 deprived23, but a dry flaky skin and an impaired epidermal barrier function— characteristics seen in AD24, may increase oxygenation and select for the low abundance of strictly anaerobic bacteria such as Lactobacillus spp or Finegoldia spp. At anaerobic conditions bacteria are fermenting organic matter, e.g. in skin the amino acid serine originating from filaggrin degradation, forming in particular lactic acid, propionic acid and other short chain fatty acids (SCFA). These metabolites lower the skin pH to pH < 5.5, maintaining a protective acidity of the skin. Potassium lactate is also one of the most relevant “natural moisturizing factors”of healthy skin. Furthermore, gram-positive anaerobe cocci such as Finegoldia,Anaerococcus, and Peptoniphilus, stimulate rapid induction of antimicrobial peptides response in human keratinocytes, which could be an important signaling mechanism to the keratinocytes when the skin is injured25. In the complete or partial absence of these organisms, danger signaling in keratinocytes and other barrier functions could be impaired, potentially favoring colonization by S. aureus. We identified AD and PSO transcriptomes which overlap substantially with previously published studies26–29. While Fig. 4 The “S. aureus signature”and functional associations. aUsing AD-associated genes (AD patients, n=82) identified at FDR level 10−5, we created an AD gene co-expression network, which was partitioned into modules by network community detection. Top enriched GO terms are indicated for each network module. bDifferential analysis between S. aureus “high”(n=27) and “low”(n=25) samples revealed 256 differentially expressed genes (FDR < 0.05, FCH ≥1.5). Hypergeometric tests revealed significant enrichment in S. aureus-associated genes (colored red) that mapped to modules M1 and M5. Gene annotations and Ingenuity Pathway analysis identified cenriched gene ontology terms in the S. aureus signature, and dpredicted upstream regulators, respectively. eMolecular networks generated between top functions and associated genes. The red color of the gene symbols indicates upregulated genes, green indicates downregulated genes. Red colored edges indicate predicted activation, yellow edges indicate inconsistent findings, and gray edges lack a predicted effect. fStatistical analysis (unpaired ttest) of RNA expression levels of selected genes in S. aureus “high”and “low”abundance samples. The center line in the dot plots corresponds to the mean, and the error bars correspond to the standard deviation. The source data files used to generate the present figure are available from the NCBI Sequence Read Archive under accession PRJNA554499, and from EBI ArrayExpress under accession E-MTAB8149 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-12253-y 8NATURE COMMUNICATIONS | (2019) 10:4703 | https://doi.org/10.1038/s41467-019-12253-y | www.nature.com/naturecommunications transcriptomic changes in PSO were dramatic, microbial changes were rather small. Yet, the opposite applied to AD lesions, suggesting a nonlinear relationship between the skin microbiota and the host transcriptome. Instead, the communication between host and microbiota might depend on defined subsets of genes, and the interaction is likely bi-directional, with the microbiota influencing host gene expression, which in turn forms a habitat for specific microbes. Previous studies have identified colonization by S. aureus in AD; however, little is known about its potential mechanistic impact10,13. We identified differentially regulated host genes between AD lesional skin samples containing high or low abundance of S. aureus, functionally enriched for mainly three activities: skin barrier function, immune activation and trp metabolism. The S. aureus high associated microbiome was overrepresented by bacterial toxins (alpha-toxin and delta-toxin), the two-component system and glucose metabolism, essentially contributed by S. aureus specific OTUs. The highly adaptable and potent pathogen S. aureus is known to express a multitude of virulence factors, including toxins that impact on the skin barrier and the immune system. The pore forming alpha-toxin, for instance, which was highly enriched in the S. aureus “high” samples, likely plays a key role in the disruption of the skin barrier in AD patients, and through being able to activate the inflammasome, resulting in the secretion of IL-1beta, alpha-toxin is able to promote further inflammation30,31. Delta-toxin, in turn, is known to contribute to AD-associated pathology through the induction of mast cell degranulation and Th2 differentiation13. The host responded to S. aureus through upregulating betadefensins (eg. DEFB4), and the expression of S100 protein family members (S100A8,S100A9,S100A7). Beta-defensins are potent antimicrobial peptides, contributing directly to local immune responses as chemoattractants for leukocytes, and through activating antigen-presenting cells (APCs)32. Psoriasin (S100A7) preferentially kills Escherichia coli, but has also bactericidal activity against S. aureus33. Being upregulated by proinflammatory cytokines, psoriasin functions as a T cell and neutrophil chemotactic agent. S100A8 and binding partner S100A9 are, in turn, expressed and released by activated phagocytes, and have powerful antimicrobial activities through the sequestration of essential trace elements. The S100A8/A9 complex also protects the host from infection by triggering toll-like receptor 4 (TLR4) and receptor for advanced glycation end-products (RAGE)- mediated inflammatory pathways, and through the recruitment of neutrophils34,35. Thus, besides targeting microorganisms, the antimicrobial defense response promotes inflammation, possibly feeding additionally into AD pathology. Further, samples abundantly colonized by S. aureus, displayed significant changes in the expression of barrier genes, such as CLDN8, a component of tight junctions (TJs). TJs are key to epithelial barriers, reside immediately under the stratum corneum, and regulate the passage of water, ions and solutes through the skin. TJ defects, including dysregulated CLDN8, are associated with AD36,37. Moreover, the expression of metalloproteinases ADAM12 and MMP12, which are implicated in tissue remodeling, cytokine and growth factor shedding, cell migration and adhesion38, was modified. Alpha-toxin is known to specifically interact ADAM10, leading to the cleavage of cadherins in epithelial cell tight-junctions, and resulting in the destruction of cell −cell contacts39, and consequently, the expression of ADAM10 was modified in S. aureus high samples. Among proinflammatory factors that were influenced by S. aureus, we note dysregulation of IL1B, CCL2, CCL19 and members of the Th2 signaling pathway. IL1B is a central mediator of inflammation40, affecting all cells of the innate immune system and playing a key role in the differentiation and function of adaptive lymphoid cells. The chemokine CCL2, which is expressed mainly in basal keratinocytes, acts through attracting monocytes, dendritic cells, Th1 and Th2 cells41, and the chemokine CCL19 is involved in the colocalization of CCR7-positive dendritic cells with CCR7-positive T cells42, facilitating antigen presentation and activation of T cells. Among Th2 associated genes, we observed the induction of chemokine receptor CCR4, cytokines IL5 and IL13, cytokine receptor IL4R, and TNFRSR4, all central mediators of allergic pathomechanisms. Finally, the AD-associated peptidase inhibitor 3 PI3/elafin, a product of Th17 activation and implicated in Th2 differentiation43, was significantly modulated in S. aureus high samples. Colonization of the skin by S. aureus activates defense and tissue repair responses, which cause increased metabolic demands, requiring the tissue to switch to glycolysis44. When densely present, S. aureus may generate localized hypoxia and ATOPIC DERMATITIS Metabolic shift to glycolysis Dominance of single microbe Shifting microbial communitees Burkholderia sp. Pelomas sp. DOWN Lactobacilli sp. Methylophilaceae g. sp S. aureus Leptotrichia sp. Tryptophan metabolism Immune activation Leukocyte migration Lactobacillus sp. HEALTHY L. crispatus ANAEROBIC NICHE PSORIASIS C. kroppenstedtti Peptostreptcoccus sp. UP C. simulans Neisseriaceae sp. Actinomyces sp. Rothia mucilaginosa Th1 pathway IFN signaling DERMIS EPIDERMIS Skin development and barrier function Streptcoccus sp. Corynebacterium sp. Pore-forming toxins Metabolic reprogramming Fig. 5 Host−microbe interaction in atopic dermatitis and psoriasis. AD is characterized by overgrowth of S. aureus, and loss of microbial diversity. In AD, colonization of the skin by S. aureus is associated with dysregulation of genes involved in epithelial barrier function, immune activation, leukocyte migration, trp degradation and metabolic reprogramming. PSO is associated with multiple species, including increased colonization by C. simulans and C. kroppenstedtii, and a loss of Lactobacillus,P. acnes and Corynebacterium spp., which may play regulatory roles NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-12253-y ARTICLE NATURE COMMUNICATIONS | (2019) 10:4703 | https://doi.org/10.1038/s41467-019-12253-y | www.nature.com/naturecommunications 9