Skin microbiome differentiates into distinct cutotypes with unique metabolic functions upon exposure to polycyclic aromatic hydrocarbons

Abstract

Background

The effects of air pollutants, particularly polycyclic aromatic hydrocarbons (PAHs), on the skin microbiome remain poorly understood. Thus, to better understand the interplay between air pollutants, microbiomes, and skin conditions, we applied metagenomics and metabolomics to analyze the effects of PAHs in air pollution on the skin microbiomes of over 120 subjects residing in two cities in China with different levels of air pollution.

Results

The skin microbiomes differentiated into two cutotypes (termed 1 and 2) with distinct taxonomic, functional, resistome, and metabolite compositions as well as skin phenotypes that transcended geography and host factors. High PAH exposure was linked to dry skin and cutotype 2, which was enriched with species with potential biodegradation functions and had reduced correlation network structure integrity. The positive correlations identified between dominant taxa, key functional genes, and metabolites in the arginine biosynthesis pathway in cutotype 1 suggest that arginine from bacteria contributes to the synthesis of filaggrin-derived natural moisturizing factors (NMFs), which provide hydration for the skin, and could explain the normal skin phenotype observed. In contrast, no correlation with the arginine biosynthesis pathway was observed in cutotype 2, which indicates the limited hydration functions of NMFs and explains the observed dry skin phenotype. In addition to dryness, skin associated with cutotype 2 appeared prone to other adverse conditions such as inflammation.

Conclusions

This study revealed the roles of PAHs in driving skin microbiome differentiation into cutotypes that vary extensively in taxonomy and metabolic functions and may subsequently lead to variations in skin–microbe interactions that affect host skin health. An improved understanding of the roles of microbiomes on skin exposed to air pollutants can aid the development of strategies that harness microbes to prevent undesirable skin conditions.

Video Abstract

Introduction

The human skin hosts a diverse community of microbes (i.e., microbiome) on the stratum corneum layer of the epidermis. The microbiome is indispensable for the modulation of host immunity, protection against opportunistic pathogens, and maintenance of skin physiology [1]. More importantly, skin conditions such as acne [2, 3] and atopic dermatitis [4, 5] have microbial contributions as part of their etiology, and skin microbiome dysbiosis may partly explain the onset of these conditions. Given the interconnectedness between the skin microbiome and host skin health, examining how host physiological and environmental factors shape the skin microbiome may allow the metabolic functions of microbes to be harnessed to protect hosts against various skin conditions. The understanding of the skin microbiome has increased in recent years owing to metagenomics sequencing [6,7,8] and mass-spectrometry-based profiling of microbial metabolites [9,10,11].

Studies have identified distinguishable microbial sub-community clusters termed cutotypes [12] (or cutaneotypes [13]) in host skin, which differ in microbial diversity, composition, functional, and antimicrobial resistance profiles. The distribution of cutotypes may vary with skin characteristics of individual hosts and thus help explain how alterations in skin microbial taxonomic composition can lead to the onset of skin conditions across different demographic groups. However, little information is currently available on factors that explain the differentiation of cutotypes and how this relates to skin health. Such knowledge would undoubtedly aid the identification of cutotype markers associated with undesirable skin conditions, thereby assisting the development of personalized solutions for improving skin health.

Chronic exposure to air pollutants is a serious concern for individuals living in urban environments. Many pollutants occur in the atmosphere, of which polycyclic aromatic hydrocarbons (PAHs) are among the most hazardous class of organic molecules to human health [14]. Earlier reports on the impacts of air pollution on skin have pointed toward an association with premature skin aging [15,16,17,18]. Specifically, cutaneous exposure to PAHs has been linked to changes in the physiological properties of the skin, potentially leading to pigmented spots [11] and cancer [19]. While specific members of the skin microbiome have been shown to metabolize PAHs [20, 21], the roles of PAH exposure in shaping the skin microbiome warrant further investigation. In our previous study [22], we demonstrated a dose–response relationship between the levels of exposure to various PAHs and changes in the taxonomic compositions of skin microbiomes, as well as variations in the abundances of microbial functional genes crucial for host–microbe interactions, virulence, and host immune modulation. However, these findings were largely based on amplicon sequencing of the phylogenetic marker genes of bacteria and fungi, while functional profiling by metagenomics sequencing was performed for only a small subset of samples with no corresponding metabolite characterization. More importantly, studies have yet to answer whether and how PAH exposure drives skin microbiome differentiation, what is the relationship between PAH exposure and cutotype differentiation, and whether any changes in cutotypes associated with PAH exposure can also result in changes in the function and metabolites of the microbiome as well as the skin–microbe interactions that may be detrimental to host skin health.

To address these questions, we performed a large-scale study of the facial cheek microbiomes of 124 Chinese females from two cities with different PAH exposure levels. The concurrent assessment of the metagenomes and metabolites allowed comprehensive characterization of the dynamic interplays between PAH exposure, skin microbiomes, and their functional potentials and metabolites. We report, for the first time, a link between PAH exposure from air pollution and the differentiation of the skin microbiome into cutotypes with distinct taxonomic, functional, resistome (i.e., the repertoire of antibiotic resistance genes), and metabolite compositions as well as skin phenotypes.

Results

Differentiation of skin microbiomes into distinct cutotypes was associated with PAH exposure

We detected 1525 species-level taxa at an average relative abundance of > 0.00001% of all the skin microbiomes, with 42 of these taxa considered core (i.e., present across all samples) as they comprised an average of nearly 80% of the reads per sample (Supplementary Fig. 1a). The microbial diversity was higher in subjects from the more polluted city of Baoding (Supplementary Fig. 1b). Fungal reads composed an average of 2.2% (95% confidence interval 0–8.7%) of all the skin microbiomes, most of which were classified as Malassezia restricta and Malassezia globosa (Supplementary Fig. 1c-d). Viral and archaeal reads cumulatively composed an average of 0.05% (95% confidence interval 0–0.36%) of the microbiome across each sample.

Prediction strength analysis estimated the presence of two distinct clusters among the samples. These two distinct community clusters, which explained 32.3% of the total community variation (Supplementary Table 1), were also detected with the Bray–Curtis dissimilarity (Fig. 1a). The microbiome of one cluster (hereafter, cutotype 1), which had a relatively lower Shannon diversity, was enriched with Cutibacterium acnes and its taxonomic relative Cutibacterium granulosum (Fig. 1a–c). Alternatively, the second cluster (hereafter, cutotype 2) had a relatively higher microbial diversity with > 100 species-level taxa (e.g., Xanthomonas citri and Rhodococcus opacus), including common skin colonizers such as Staphylococcus aureus and Corynebacterium species (Fig. 1a–c).

Fig. 1

Cheek microbiomes were differentiated into cutotypes and associated with PAH exposure. a Principal coordinate analysis (PCoA) plot based on the Bray–Curtis dissimilarity of microbial community composition. Data points represent metagenomic samples with point size representing Shannon diversity and color showing the cutotype as determined by the prediction strength analysis. The inset panel shows the pairwise Bray–Curtis dissimilarity between samples of the same or different cutotypes, showing that microbiomes within a cutotype were more similar than between cutotypes (**Tukey’s pairwise post-hoc p < 0.05 for all pairwise comparisons). b PCoA plot as shown in panel a with the data points colored according to the relative abundance of Cutibacterium acnes, Xanthomonas citri, and Rhodococcus opacus in cutotypes 1 and 2. c Top 12 species-level bacterial taxa based on the average relative abundance in all the samples organized by cutotype. Other bacterial taxa were grouped into the “Minor/Unclassified” category. d Distance-based redundancy analysis depicting the effects of exposure to high and low molecular weight PAHs (purple and green, respectively) on cutotype differentiation. PAHs with a molecular weight < 200 g/mol were considered low. e Spearman’s correlation between the first component of the PCoA and the measured concentrations of the six PAHs. Linear regression lines and R² values are shown. f Association of samples in the respective cutotypes with different host factors and host skin phenotype

Cutotype 2 was enriched with a diverse array of genera (e.g., Paracoccus, Archromobacter, Caulobacter, Pseudomonas, Rhodococcus, and Sphingomonas) known to have biodegradation capabilities (Supplementary Table 2) and some of them were found in PAH-degrading microbial communities [23, 24]. We thus hypothesized that the differentiation of the microbiome into cutotype 2 was associated with PAH exposure. Indeed, distance-based redundancy analysis revealed that exposure to many different PAHs appeared to drive the microbiome toward cutotype 2 (Fig. 1d). Exposure to a small number of lower-molecular-weight PAHs (i.e., < 200 g/mol) was associated with the differentiation of all cutotype 1 microbiomes and a small subset of cutotype 2 microbiomes, whereas exposure to a large number of higher-molecular-weight PAHs (i.e., > 200 g/mol) was associated with the differentiation of most cutotype 2 microbiomes. The exposure levels to many PAHs were also correlated with the first dimension of the community compositional principal coordinate axis (i.e., PCoA1 in Fig. 1a) (Fig. 1e), suggesting that PAH exposure and community composition are partially linked.

Cutotype was not related to age group and only marginally associated with city (Fig. 1f, Supplementary Fig. 1e and Supplementary Table 1). Subjects with dry skin (sebum level < 70 μg/cm²) were prone to having a microbiome resembling cutotype 2, while those with normal to high sebum levels (> 70 µg/cm²) were prone to having a microbiome resembling cutotype 1 (Fig. 1f, chi-square with Yates correction = 14.4, p = 0.00015), which was consistent with the high abundance of Cutibacterium. Although not statistically significant, subjects with higher pigmentation frequency were more prone to having a microbiome resembling cutotype 2 (Fig. 1f). Regarding fungal taxa, community differences could neither be explained by cutotypes nor other host factors, and shifts along PCoA1 based on the Bray–Curtis dissimilarity were only associated with exposure to nicotine (Spearman’s ρ = 0.308, p = 0.018) and cotinine (Spearman’s ρ = 0.251, p = 0.045).

Species-specific correlations between growth rate and PAH exposure

Growth rate inference [6, 7] was performed to assess how cutotypes and PAH exposure may influence microbial physiology. Our findings suggested that taxa in the genera of Micrococcus and Corynebacterium were among those with the highest inferred growth rate and it was similar between the two cutotypes (Fig. 2a). Taxa within Enhydrobacter, a genus phylogenetically related to Moraxella and hypothesized to be enriched on Asian skin [12, 25, 26], were also inferred to have a high growth rate. The levels of exposure to dibenzo[a,h]anthracene and acenaphthene were associated with changes in the inferred growth rates of two bacterial species (i.e., a bacterium affiliated with Micrococcaceae and Janibacter indicus) (Supplementary Table 2). City, acne onset, and pigmentation frequency were not associated with changes in the inferred growth rates of any bacterial taxa. Similarly, single-nucleotide polymorphism-based subspecies-level growth rate estimation [9] did not indicate cutotype-based growth pattern differences between the clusters of related taxa within C. acnes (Fig. 2b) and Micrococcus luteus (Fig. 2c). Inferred growth rates of C. acnes and M. luteus subspecies taxa on the skin were comparable to those reported in previous studies [6, 7, 27]. Overall, pollution exposure did not appear to influence inferred growth rates across the microbiome.

Fig. 2

Growth rate inference of the cheek microbiome at the species and subspecies levels. a Heatmap of the species level inferred growth rates for taxa across samples grouped by cutotype. b, c Heatmaps of the inferred growth rates for subspecies taxa belonging to b Cutibacterium acnes and c Micrococcus luteus

PAH-associated cutotypes were enriched in diverse potential functions

The functional profile was examined to understand how compositional differences between cutotypes translated to variations in functional potentials. The congruence between taxonomy and functional pathways in the skin microbiomes in cutotypes 1 and 2 (Fig. 3a) suggested that samples with similar taxonomic compositions also tended to have similar functions. As in taxonomy, cutotype was the strongest factor in explaining functional differences between microbiomes, as compared to other parameters such as clinical or city factors (Supplementary Table 1).

Fig. 3

Cutotypes differed in contributional functional diversity and enriched functional pathways. a Procrustes analysis showing congruence between the taxonomic and functional compositions of the microbiome. The statistically significant Pearson correlation indicates that microbiomes were similar in both taxonomic composition and functions. b Contribution diversity differences between cutotypes of core gene UniRef90 pathways (those present in ≥ 75% of the samples within each cutotype). Core pathways detected in both cutotypes are considered to be “shared” (triangles). c Examples of core pathways (folate transformation II and L-valine biosynthesis) show differences in contributional diversity between cutotypes. Relative abundance is expressed in reads per kilobase per million reads (RPKM). d Enrichment of UniRef90 pathways between cutotypes according to MaAslin2. Only statistically significant (adjusted p < 0.05) pathways are shown. e Relative abundance in RPKM of KOs involved in xenobiotic degradation across samples in the respective cutotypes

The core set of UniRef90 pathways (i.e., ≥ 75% of samples in each cutotype) was assessed by contributional diversity analysis to further compare the pathway diversity between cutotypes. For both cutotypes, the core was dominated by pathways involved in the biosynthesis and metabolism of nucleotides and amino acids (Supplementary Table 3). Despite this conservation, core pathways for cutotype 1 mainly exhibited low within- and between-sample diversity (i.e., “simple and conserved” contributional diversity [28]) (Fig. 3b). In contrast, core pathways in cutotype 2 were dominated by those with high between sample diversity (i.e., “simple and variable” and “complex and variable” contributional diversity [28]) (Fig. 3b). The discrepancy between the two cutotypes regarding within and between sample variations was most likely due to the dominant contribution of C. acnes in cutotype 1 to the shared core pathways (e.g., folate transformation and valine biosynthesis) (Fig. 3c).

Multivariate analysis was performed to detect the enrichment of KEGG Orthology (KO) gene families between cutotypes, which revealed the enrichment of 556 KOs (Supplementary Table 2) and 53 UniRef pathways spanning multiple facets of microbial physiology (Fig. 3d) in cutotype 2 versus cutotype 1. At the pathway level, cutotype 1 was characterized by amino acid degradation potential (PWY_5028: L-histidine degradation III), while cutotype 2 was associated with multiple amino acid biosynthesis pathways (n = 8), including tryptophane (TRPSYN_PWY: L-tryptophan biosynthesis) and methionine biosynthesis (n = 4). Cutotype 2 was also characterized by an increased number of pathways related to carbohydrate biosynthesis (n = 6) and the tricarboxylic acid (TCA) cycle (n = 7), suggesting a shift in the microbial metabolism.

Given that high PAH exposure was linked to microbiome differentiation into cutotype 2 concomitant with the enrichment of pathways related to diverse functions, we hypothesized that this cutotype presented an increased ability to degrade xenobiotic compounds. To test this hypothesis, KOs associated with xenobiotic compound degradation were compared between the two cutotypes. The Mann–Whitney test results showed an increased abundance of KOs involved in the degradation of benzoate (FDR-adjusted p = 6.6 × 10^–16), chloroalkane/alkene (FDR-adjusted p = 7.2 × 10^–16), PAH (FDR-adjusted p = 6.9 × 10^–5), steroids (FDR-adjusted p = 9.9 × 10^–5), styrene (FDR-adjusted p < 7.8 × 10^–10), and toluene (FDR-adjusted p < 5.7 × 10^–14) in cutotype 2 (Fig. 3e). Conversely, abundances of KOs linked to nitrotoluene degradation were increased in cutotype 1 (FDR-adjusted p = 3.78 × 10^–5). Overall, these results suggested that while each cutotype was characterized by distinct microbes that contributed to shared core functions, they also presented differential enrichment of metabolic functions that conferred adaptive and biodegradative advantages given certain environmental conditions.

Cutotypes also had different resistomes, which may explain cutotype-specific responses to host antimicrobial treatment. The most abundant gene family markers were those conferring resistance to fluoroquinolones, aminoglycosides, and beta-lactams (Supplementary Table 4). Cutotype explained nearly 18% of resistome variation (Supplementary Table 1), with point mutations in gyrA (ARO: 3,003,974) conferring resistance to fluoroquinolone and an efflux protein (tet(V)) conferring resistance to tetracycline (ARO: 3,000,181) enriched in cutotype 2. Exposure to PAHs also positively correlated with the abundance of a small number of antibiotic resistance gene families (Supplementary Table 2), including those encoding for a tetracycline efflux pump (tet(33)) and several plasmid- or integron-encoded gene families related to aminoglycoside resistance. Thus, these results suggest that chronic exposure to PAH may favor the presence of species with antimicrobial resistance genes.

Skin metabolites were specific to cutotypes and correlated with the microbiome and PAHs

In our earlier report using 16S rRNA gene amplicon sequencing data [11], we showed that the skin metabolome accounted for one third of the variability in bacterial diversity. Here, we further investigated the links between metabolome and metagenomic data. The metabolomic profiles of subjects [11] were examined with the corresponding metagenomes to test whether the molecules on skin were specific to cutotypes. Many gamma-glutamyl amino acids, xenobiotics, and unknown metabolites were overrepresented in cutotype 2 (Supplementary Fig. 2). In contrast, many diverse lipid molecules, largely represented by the sebum components myristic and palmitoleic acid [29], were overrepresented in cutotype 1.

Correlation network analysis was performed by integrating metabolomic data with corresponding taxonomic, functional, metabolomic, and PAH features for each cutotype to (i) provide a comprehensive and multifaceted view of the potential ecological relationships between pollutant exposure and microbial and metabolic features and (ii) compare the relationships between microbiome components of the two cutotypes. While networks of the cutotypes presented similar numbers of features (i.e., nodes), the network of cutotype 1 was more integrated with a higher number of edges and average number of neighbors per node (Supplementary Table 5). The network density of cutotype 1 also exceeded that of cutotype 2, which suggested a more stable network in the former.

To observe the overall correlations between the microbiome features of the two cutotypes, their networks were merged, which produced 2684 nodes and 52,727 correlations (Supplementary Table 5). Significant correlations involving metabolites previously [11] shown to be associated with undesirable skin phenotypes (e.g., N-acetyl and G-glutamyl amino acids, urea cycle intermediates, amino acids of tryptophan metabolism, fatty acids and carnitine lipids, lactate, and TCA cycle metabolites) were observed (Fig. 4). N-acetylglutamine, a metabolite of arginine biosynthesis, was central to the network and had significant positive correlations with multiple Corynebacterium species (n = 43) and KOs (n = 30) (e.g., K01607 that encodes 4-carboxymuconolactone decarboxylase, an enzyme involved in benzoate degradation) in cutotype 2 (Supplementary Table 5 and Fig. 4). N-acetylglutamine was also correlated with the fpr gene (K00528), which encodes an oxidoreductase, as well as KOs of the TCA cycle, carbohydrate metabolism, and biofilm formation (Supplementary Table 5). These functions are representative of the bacterial degradation strategy for aromatic compounds, in which compounds are de-aromatized and then funneled into central carbon pathways such as the TCA cycle [23]. The strong positive correlations with Corynebacterium species suggested that the cluster may contribute to PAH biodegradation. Other species capable of biodegradation, including Pseudomonas and Rhodococcus, were also positively correlated with N-acetylglutamine (Supplementary Table 5). For cutotype 1, two clusters were observed. One cluster featured species, including C. acnes, and KOs that were positively correlated with N-delta-acetylornithine, which is another arginine pathway metabolite. The second cluster showed positive correlations with 2-oxoarginine and many KOs (n = 56) encoding base and mismatch repair and related functions (K03469, K01246, K03601, and K01971). 2-Oxoarginine was also positively correlated with glutaryl-CoA dehydrogenase (K00252) for benzoate degradation and glutamate dehydrogenase (K00262) for arginine biosynthesis (Fig. 4).

Fig. 4

Cutotype-specific correlations between compositional, functional, metabolomic features, and PAH exposure. Each node represents a feature and is colored by feature type and sized by node degree. Blue and orange edges represent those in cutotypes 1 and 2, respectively. Solid and dashed edges represent positive and negative correlations, respectively, with line thickness representing the absolute value of the Spearman’s correlation. Only metabolites considered enriched in the polluted city of Baoding [11] are shown to enhance visual clarity. Nodes connected to N-acetylglutamine are circled to highlight a cluster of Corynebacterium taxa. The species-level taxonomy of these Corynebacterium nodes, as well as all other significant correlations, can be found in Supplementary Table 5. The multivariate analysis was performed using HallA. Correlations with corrected q value ≤ 0.25 are statistically significant [30], and only those with a Spearman’s correlation >|0.25| were included in analysis

Kynurenate, a tryptophan pathway metabolite, which may stimulate the onset of inflammatory responses in hosts [31], showed positive correlations with the oral bacterium Streptococcus sobrinus in cutotype 2. Interestingly, the presence of oral bacteria on the skin has been attributed to skin aging [32]. An unknown metabolite, X-13737, previously shown to be enriched in subjects from the more polluted city of Baoding [11], showed a positive correlation in cutotype 2 with a taxon affiliated with Pseudomonas, which has been postulated to compete with resident skin Cutibacterium [33].

Discussion

The skin microbiome is shaped by the external environment surrounding the skin [34], and subsequently, exposure to pollutants such as PAHs has been linked to a wide range of skin conditions [14, 19, 35, 36]. Recently, the roles of the skin microbiome in modulating host responses to environmental stressors have gained recognition [37]. By analyzing 32 metagenomic samples from cheeks and scalps, we previously associated PAH exposure with changes in the compositional and functional capabilities of the skin microbiome, with potential consequences for skin phenotype [22]. In this study, by analyzing the metagenomes and corresponding metabolites from the cheeks of over 120 subjects, we provided evidence that PAH exposure could play a role in the differentiation of skin microbiomes into cutotypes that have distinct compositions and metabolic functions at the community and genomic levels and are associated with different skin phenotypes. The link between PAH exposure and cutotype differences may help explain the potential roles of the microbiome in modulating host skin health upon air pollutant exposure.

The presence of cutotypes on the skin of subjects in this study was consistent with the stratification of microbiomes on other human body parts, which has been linked to changes in host phenotypes, physiology, and susceptibility to diseases [38,39,40]. A recent study of nearly 300 Chinese individuals also reported that skin microbiomes can be grouped into two cutotypes [12], with the differentiation mainly driven by nutrient availability as well as host physiology and chronological age [8]. Our study showed that cutotype differentiation may also occur as a function of environmental exposure to air pollution. The bacterial taxa previously documented to have biodegradation capabilities (e.g., Archromobacter, Caulobacter, Paracoccus, Pseudomonas, Rhodococcus, and Sphingomonas) [23, 24] were enriched in cutotype 2, and a positive correlation was detected between acenaphthene exposure level and the inferred growth rate of a species in Janibacter, a genus capable of using multiple PAH types as carbon sources [41]. These observations together with the documented PAH metabolic capabilities of common skin colonizers [20, 21] suggest that PAHs are actively metabolized by the skin microbiome. The species with biodegradation abilities may have been originated from the environment [24, 42] and managed to transiently colonize on skin due to the presence of PAHs. However, PAHs do not appear to greatly influence fungal community composition. The cutotypes also showed significant resistome differences, which was consistent with previous cutotype analyses [12]. The roles of PAH exposure in driving resistome variations via cutotype differentiation remain to be elucidated.

Multivariate analysis of metagenomics, metabolomics, and PAH exposure data revealed complex associations between pollutant exposure and various aspects of the skin microbiome. The correlation network of cutotype 2 was strongly associated with exposure to multiple PAHs and had a reduced integrity, which has been previously associated with acne and dandruff on the cheeks and scalp [22, 43]. However, whether microbiome network stability can be used as a diagnostic marker to test for susceptibility to pollutant-related skin conditions requires further investigation. Each PAH appeared to have unique and cutotype-specific associations with microbial features. Metabolites previously detected at higher levels in the more polluted city [11] showed cutotype-specific associations with functional genes involved in DNA repair, secretion systems, and pollutant degradation. While we cannot determine whether the detected metabolites on the skin originated from microbial metabolism, environmental exposure, or anthropogenic factors [10], the presence of shared but opposing correlations between cutotypes suggest that specific microbes interact differently with chemicals on the skin.

The two identified cutotypes were associated with distinct skin dryness and microbial metabolism in which the arginine pathway appeared to play a pivotal role (Fig. 5). Cutotype 1, which was mostly associated with normal to greasy skin, was dominated by C. acnes (consistent with previous reports on sebaceous sites [8]) and functions related to lipid metabolism and amino acid biosynthesis. C. acnes, among other taxa, was positively correlated with N-delta-acetylornithine as well as many KOs (i.e., K00619 [argA], K00930 [argB], K00145 [argC], K00611 [argF], K01755 [argH], and K00620 [argJ]) in the arginine biosynthesis pathway, suggesting that arginine is produced in cutotype 1. Arginine is a key amino acid in filaggrin-derived natural moisturizing factors (NMFs), which provide hydration for the skin [44]. We thus speculate that arginine from bacteria contributes to skin hydration (Fig. 5). High-throughput mass spectrometry proteomic analysis of the cheek stratum corneum of a subset of subjects (34 from cutotype 1 and 19 from cutotype 2) revealed that the expression level of the FLG gene-encoded filaggrin [45] was similar between the two cutotypes (data not shown), suggesting that skin dryness does not appear to be linked to the FLG gene.

Fig. 5

Schematic illustration of the potential influence of the arginine biosynthesis pathway in cutotypes 1 and 2 on natural moisturizing factors (NMFs) for skin hydration. Arginine biosynthesis played an important role in cutotype 1, but only a limited role in cutotype 2. The illustration was drawn based on the key correlations found between the dominant taxa, metabolic functions, and metabolites in the respective cutotypes. The KOs and metabolites highlighted by a rectangular box are those identified as having significant positive correlations in the network analysis

In contrast, cutotype 2 was associated with dry skin and hyperpigmentation in subjects aged under 45. Diverse species, including Paracoccus [46], Sphingomonas [23], and Pseudomonas [20], that can reportedly degrade PAHs were found in cutotype 2. In particular, members of Corynebacterium were the key taxa in cutotype 2. These taxa were positively correlated with benzoate degradation and KOs in the TCA cycle (i.e., K01607 and K00262 [gdhA]). Growth in the presence of aromatic compounds has been shown to upregulate enzymes in the TCA cycle of Corynebacterium glutamicum (e.g., a fivefold increase in GdhA that leads to L-glutamate production) [47], which we speculate is similar to that in the microbiome of skin exposed to PAHs. Consistent with the potential increase in L-glutamate, the metabolite N-acetylglutamine, produced from L-glutamate, was positively correlated with Corynebacterium taxa (Fig. 5). In contrast to cutotype 1, Corynebacterium taxa in cutotype 2 were not correlated with any KOs or metabolites in the arginine biosynthesis pathway, and positive correlations with KOs in the arginine biosynthesis pathway were only found in taxa with low relative abundances. Thus, arginine synthesis by the microbiome may be limited in cutotype 2 compared with cutotype 1, which indicates the limited hydration functions of NMFs and explains the dry skin phenotype associated with cutotype 2 (Fig. 5). In addition to dryness, skin associated with cutotype 2 may be prone to other adverse conditions because taxa and KOs in cutotype 2 were also positively correlated with kynurenate, a key tryptophan derivative linked to pathological processes such as atopic dermatitis, which causes dry, itchy, and inflamed skin [48, 49]. Kynurenate could be derived from microbes to trigger host-microbe interactions as it is a ligand for the transcription factor aryl hydrocarbon receptor that regulates the functions of many skin cell types [50]. With increasing recognition of hydrocarbons in polluted air as an important factor in the development of atopic dermatitis, particularly in childhood [51,52,53], our results suggest that shifts in skin microbiome composition after air pollution exposure could contribute to the aggravation of skin dryness and inflammation.

While this study improved our understanding of the skin microbiome in the context of PAH exposure, it has some limitations. First, although the growth rates of microbial taxa can be inferred computationally, whether the taxa are alive or dead cannot be discerned from genomic DNA data alone. Second, there is little information on whether detected metabolites were derived from microbes, the host, or the environment. Third, multivariate analysis only reveals correlations but not causations such that observed associations may be ecologically significant or arise from stochastic processes [54]. Subsequently, the causal mechanisms underlying detected associations (e.g., adverse skin conditions) requires future experimental verification, such as by prospecting specific microbial–skin interplays using in vitro skin models colonized with microbial species from each cutotype and exposing them to different types and concentrations of PAHs [55,56,57]. Furthermore, multi-omics analyses including metatranscriptomics and metaproteomics together with metagenomics and metabolomics would provide finer understanding of the relationships between PAH exposure and microbial physiology, and the combination of culture-independent and culture-dependent assays would help to validate in silico associations [58].

Conclusions

This study revealed the associations between PAH exposure and differentiation of skin microbiome into distinct cutotypes. Some of the differences in cutotypes may result in changes in microbial adaptation and skin–microbe interactions that affect host skin health. The coupling of PAH exposure and cutotype differentiation will help future studies to compartmentalize skin microbiomes into clusters. In-depth characterization of these clusters can then be performed to identify taxonomic, functional, or metabolic biomarkers associated with changes in skin phenotypes (e.g., acne onset, pigmentation frequency, and aging). Such a workflow will help develop strategies for preventing pollutant-associated skin effects by identifying specific features and interactions from microbiome clusters linked to pollutant exposure.

Materials and methods

Subject characteristics

The analyzed cheek samples originated from a cohort of female participants aged 25–45 years from the Chinese cities of Baoding (more severe air pollution) and Dalian (less severe) as part of a multi-objective study examining PAH accumulation in hair [59] and the associations of PAH exposure with signs of facial aging [36], skin microbial taxonomic composition [22], and metabolomes [11]. The participants resided in the respective cities for at least 15 years at the time of sampling. Details regarding cohort recruitment, inclusion and exclusion criteria, and sample collection procedures have been previously described [22]. For each subject, comprehensive skin physiological data collected through self-reporting and clinical assessment were also previously reported [22].

Sequencing and quality control

Genomic DNA was extracted from sampling swabs using the PowerSoil DNA isolation kit (MO BIO Laboratories, Carlsbad, CA, USA) following the manufacturer’s instructions with minor modifications [22]. Library preparation and sequencing on an Ilumina NovaSeq platform to generate 150-bp paired-end reads were performed by SeqMatic LLC (Fremont, CA, USA) according to the manufacturer’s instructions. We processed and sequenced 124 samples (Baoding n = 61, Dalian n = 63) and four negative controls (i.e., new sterile swabs). The sequenced samples were selected across eight pollution exposure groups from weak to strong as previously described [11, 22] with 7 to 21 samples per group to obtain a balanced representation among exposure levels. Two samples from Baoding did not yield sufficient reads and were removed from downstream analysis. Sequence adapters were removed from the raw reads using AdapterRemoval (v2.3.1) [60] and quality filtering and human DNA read removal was performed using Kneaddata (v0.7.4) as previously described [22]. Sequencing quality information for the samples is shown in Supplementary Table 6. Following taxonomic classification of the reads using Kraken2 (v2.0.7-beta; kmer length = 35, confidence score threshold = 0) [61], then species-level abundance estimation using Bracken (v2.6.1; threshold for filter = 0) [62], putative contaminants were identified with the R package “decontam” based on the prevalence mode with 0.1 as the significance threshold [63]. This method deemed 65 species to be contaminants (Supplementary Table 6) and reads that could be mapped to the representative genomes of the contaminant species (Supplementary Table 6) were removed using an in-house script. Following quality filtering and contaminant removal, an average of 9,138,597 paired-end reads were retained per sample.

Diversity and community composition analyses

Taxonomic classification of the filtered and contaminant-free reads was performed using Kraken2 and Bracken as described above. The species-level taxonomy was used to identify host factors associated with community diversity and compositional changes. Based on the taxonomic data, prediction strength analysis to estimate the optimal number of clusters in the community was performed based on the “partitioning around medoids” discrete clustering method [64] using the “prediction.strength” command (with 100 random splits) in the R package “fpc” (v2.2–9). The highest mean prediction strength was 0.81 when k = 2, where k is the number of clusters. This result formed the basis for the presence of two cutotypes among all samples. FindFungi (v0.23) [65] was used to assign taxonomy to the members of the fungal communities.

HUMAnN3 (v3.0.0.alpha.3) [66] was used to profile the potential metabolic functions of the metagenomes, revealing 507,517 UniRef90 gene families that were grouped into 361 MetaCyc functional pathways and 5564 KOs. Based on this UniRef90 pathway data, contributional diversity (Gini–Simpson for within-sample and Bray–Curtis for between-sample diversity) [28] was calculated for each cutotype and compared by selecting pathways detected in > 75% of samples within each cutotype. Antimicrobial resistance gene family markers were detected using ShortBRED (v0.9.3) [67]. The associations of the taxonomic, functional, and antimicrobial resistance data with host and environmental factors, including PAH exposure, were determined using MaAsLin2 [68] with the generalized linear model, with city, acne onset, wrinkle grade, facial pigmentation frequency, and exposure to PAHs as fixed effects and BMI and age as random effects. An adjusted p value (i.e., q value) ≤ 0.05 was considered statistically significant in this analysis. The composition differences (ß-diversity) between samples were analyzed based on the Bray–Curtis dissimilarity using the function “vegdist” in the R package “vegan,” and the permutational multivariate analysis of variance (PERMANOVA) test was applied using the “adonis” function in the R package “vegan” with 999 permutations to test the influence of factors including cutotype, city, acne onset, age group, and facial pigmentation frequency. The effects of PAHs on community differentiation by cutotype were assessed with distance-based redundancy analysis using the “capscale” function in the R package “vegan.” The Procrustes test was performed to determine congruency between taxonomic (Bracken) and functional (HUMAnN) composition data using the “protest” function in the R package “vegan.” To assess α-diversity, the samples were first subsampled using seqtk (v1.3) to a read depth of 313,504 reads per sample, which corresponded to the sequencing depth of the sample with the lowest number of reads. This depth was representative of the total species-level taxon richness of the metagenomes (Supplementary Fig. 3). Shannon diversity was calculated for the subsampled dataset in “vegan.” Statistical significance of comparisons between two and more than two groups were performed using the Mann–Whitney and Kruskal–Wallis tests, respectively, in the R package “stats” (v3.6.1). False discovery rate (FDR) adjustment of the statistical significance was performed using the Benjamini–Hochberg method.

Community- and strain-level in situ growth rate estimation

Growth rate estimation of the community was performed using the Growth Rate Index (GRiD) algorithm [6] with a coverage cutoff (-c) of 0.2, and ambiguous reads were re-assigned using Pathoscope (v2.0) [69]. Growth rate estimation of subspecies was performed using the Strain-Level Metagenomic Estimation of Growth Rate (SMEG) (v1.1.5) algorithm [7]. Reference genomes for each species included in the SMEG analysis (Supplementary Table 7) were obtained from the RefSeq archive of NCBI (accessed 2 June 2020).

Representation of metabolites in cutotypes

The representation of metabolites in each cutotype was characterized using a v test [11] with the “catdes” function in the “FactoMinerR” package. For continuous variables, we tested if the mean of a particular subgroup was different from the mean of the total population. For discrete or qualitative variables, we tested if the proportion of a modality in a particular subgroup was over-expressed or under-expressed compared with the total population.

Correlation analysis of the multivariate datasets

For each cutotype, taxonomic (Bracken) and functional (HUMAnN) features of the microbiome as well as the measured skin metabolites and PAHs were incorporated into a multivariate analysis to identify the correlating features between the datasets. Residualized values of all features from each type of data were determined based on a linear mixed-effects model using the R package “lme4.” The model included city, acne onset, and facial pigmentation frequency as fixed effects and age group as a random effect. The resulting residualized datasets were subjected to pairwise hierarchical all-against-all association (HAllA) analysis (v0.8.18) as previously described [30]. Spearman’s rank correlations with a Benjamini–Hochberg FDR-corrected p value (i.e., q value) ≤ 0.25 were considered statistically significant [30]. Spearman’s correlations <|0.25| were considered weak and were removed from the output. A sub-network containing 500 correlations of the lowest q value within each HAllA analysis between two datasets was constructed for each cutotype. The sub-networks for the two cutotypes were merged and visualized in Cytoscape (v3.8.2) to identify both correlations shared between the cutotypes and cutotype-specific correlations [70].

Measurement of PAHs in hair samples

The detailed methods for analyzing PAHs in hair samples to assess subjects’ pollutant exposure levels have been previously reported [59]. Briefly, hair samples were washed to remove any external deposits and then pulverized, hydrolyzed, extracted, and analyzed using gas and liquid chromatography coupled with tandem mass spectrometry (MS/MS). We quantified 15 PAHs (all part of the US-EPA priority list), nicotine, and cotinine.

Measurement of skin metabolites

Detailed descriptions of the methods used for untargeted metabolomics analysis of skin samples have been previously described [11]. Briefly, the samples were extracted with methanol and then divided into four equal fractions. Two fractions were analyzed with separate reverse-phase (RP) ultra-performance liquid chromatography (UPLC)-MS/MS with positive ion mode electrospray ionization (ESI) optimized for hydrophilic and hydrophobic compounds, respectively. One fraction was analyzed with RP/UPLC-MS/MS with negative ion mode ESI, and one was analyzed with hydrophilic interaction liquid chromatography/UPLC-MS/MS with negative ion mode ESI. By comparing against a library of pure standards or using the analytical profiles (i.e., retention time, molecular weight, preferred adducts, and in-source fragments), 468 metabolites were identified. Those that could not be matched to known compounds were indicated with an X followed by a number.