The impact of skin care products on skin chemistry and microbiome dynamics

Skin care and hygiene products persist on the skin

Systematic strategies to influence both the skin chemistry and microbiome have not yet been investigated. The outermost layer of the skin turns over every 3 to 4 weeks [8, 9]. How the microbiome and chemistry are influenced by altering personal care and how long the chemicals of personal care products persist on the skin are essentially uncharacterized. In this study, we collected samples from skin of 12 healthy individuals—six males and six females—over 9 weeks. One female volunteer had withdrawn due to skin irritations that developed, and therefore, we describe the remaining 11 volunteers. Samples were collected from each arm, armpit, foot, and face, including both the right and left sides of the body (Fig. 1a). All participants were asked to adhere to the same daily personal care routine during the first 6 weeks of this study (Fig. 1b). The volunteers were asked to refrain from using any personal care product for weeks 1–3 except a mild body wash (Fig. 1b). During weeks 4–6, in addition to the body wash, participants were asked to apply selected commercial skin care products at specific body parts: a moisturizer on the arm, a sunscreen on the face, an antiperspirant on the armpits, and a soothing powder on the foot (Fig. 1b). To monitor adherence of participants to the study protocol, molecular features found in the antiperspirant, facial lotion, moisturizer, and foot powder were directly tracked with mass spectrometry from the skin samples. For all participants, the mass spectrometry data revealed the accumulation of specific beauty product ingredients during weeks 4–6 (Additional file 1: Figure S1A-I, Fig. 2a orange arrows). Examples of compounds that were highly abundant during T4–T6 in skin samples are avobenzone (Additional file 1: Figure S1A), dexpanthenol (Additional file 1: Figure S1B), and benzalkonium chloride (Additional file 1: Figure S1C) from the facial sunscreen; trehalose 6-phosphate (Additional file 1: Figure S1D) and glycerol stearate (Additional file 1: Figure S1E) from the moisturizer applied on arms; indolin (Additional file 1: Figure S1F) and an unannotated compound (m/z 233.9, rt 183.29 s) (Additional file 1: Figure S1G) from the foot powder; and decapropylene glycol (Additional file 1: Figure S1H) and nonapropylene glycol (Additional file 1: Figure S1I) from the antiperspirant. These results suggest that there is likely a compliance of all individuals to study requirements and even if all participants confirmed using each product every day, the amount of product applied by each individual may vary. Finally, for weeks 7–9, the participants were asked to return to their normal routine by using the same personal care products they used prior to the study. In total, excluding all blanks and personal care products themselves, we analyzed 2192 skin samples for both metabolomics and microbiome analyses.

To understand how long beauty products persist on the skin, we monitored compounds found in deodorants used by two volunteers—female 1 and female 3—before the study (T0), over the first 3 weeks (T1–T3) (Fig. 1b). During this phase, all participants used exclusively the same body wash during showering, making it easier to track ingredients of their personal care products. The data in the first 3 weeks (T1–T3) revealed that many ingredients of deodorants used on armpits (Fig. 2a) persist on the skin during this time and were still detected during the first 3 weeks or at least during the first week following the last day of use. Each of the compounds detected in the armpits of individuals exhibited its own unique half-life. For example, the polyethylene glycol (PEG)-derived compounds m/z 344.227, rt 143 s (Fig. 2b, S1J); m/z 432.279, rt 158 s (Fig. 2b, S1K); and m/z 388.253, rt 151 s (Fig. 2b, S1L) detected on armpits of volunteer 1 have a calculated half-life of 0.5 weeks (Additional file 1: Figure S1J-L, all p values < 1.81e−07), while polypropylene glycol (PPG)-derived molecules m/z 481.87, rt 501 s (Fig. 2c, S1M); m/z 560.420, rt 538 s (Fig. 2c, S1N); m/z 788.608, rt 459 s (Fig. 2d, S1O); m/z 846.650, rt 473 s (Fig. 2d, S1P); and m/z 444.338, rt 486 s (Fig. 2d, S1Q) found on armpits of volunteers 3 and 1 (Fig. 2a) have a calculated half-life ranging from 0.7 to 1.9 weeks (Additional file 1: Figure S1M-Q, all p values < 0.02), even though they originate from the same deodorant used by each individual. For some ingredients of deodorant used by volunteer 3 on time 0 (Additional file 1: Figure S1M, N), a decline was observed during the first week, then little to no traces of these ingredients were detected during weeks 4–6 (T4–T6), then finally these ingredients reappear again during the last 3 weeks of personal product use (T7–T9). This suggests that these ingredients are present exclusively in the personal deodorant used by volunteer 3 before the study. Because a similar deodorant (Additional file 1: Figure S1O-Q) and a face lotion (Additional file 1: Figure S1R) was used by volunteer 3 and volunteer 2, respectively, prior to the study, there was no decline or absence of their ingredients during weeks 4–6 (T4–T6).

Polyethylene glycol compounds (Additional file 1: Figure S1J-L) wash out faster from the skin than polypropylene glycol (Additional file 1: Figure S1M-Q)(HL ~ 0.5 weeks vs ~ 1.9 weeks) and faster than fatty acids used in lotions (HL ~ 1.2 weeks) (Additional file 1: Figure S1R), consistent with their hydrophilic (PEG) and hydrophobic properties (PPG and fatty acids) [25, 26]. This difference in hydrophobicity is also reflected in the retention time as detected by mass spectrometry. Following the linear decrease of two PPG compounds from T0 to T1, they accumulated noticeably during weeks 2 and 3 (Additional file 1: Figure S1M, N). This accumulation might be due to other sources of PPG such as the body wash used during this period or the clothes worn by person 3. Although PPG compounds were not listed in the ingredient list of the shampoo, we manually inspected the LC-MS data collected from this product and confirmed the absence of PPG compounds in the shampoo. The data suggest that this trend is characteristic of accumulation of PPG from additional sources. These could be clothes, beds, or sheets, in agreement with the observation of these molecules found in human habitats [27] but also in the public GNPS mass spectrometry dataset MSV000079274 that investigated the chemicals from dust collected from 1053 mattresses of children.

Temporal molecular and bacterial diversity in response to personal care use

To assess the effect of discontinuing and resuming the use of skin care products on molecular and microbiota dynamics, we first evaluated their temporal diversity. Skin sites varied markedly in their initial level (T0) of molecular and bacterial diversity, with higher molecular diversity at all sites for female participants compared to males (Fig. 3a, b, Wilcoxon rank-sum-WR test, p values ranging from 0.01 to 0.0001, from foot to arm) and higher bacterial diversity in face (WR test, p = 0.0009) and armpits (WR test, p = 0.002) for females (Fig. 3c, d). Temporal diversity was similar across the right and left sides of each body site of all individuals (WR test, molecular diversity: all p values > 0.05; bacterial diversity: all p values > 0.20). The data show that refraining from using beauty products (T1–T3) leads to a significant decrease in molecular diversity at all sites (Fig. 3a, b, WR test, face: p = 8.29e−07, arm: p = 7.08e−09, armpit: p = 1.13e−05, foot: p = 0.002) and bacterial diversity mainly in armpits (WR test, p = 0.03) and feet (WR test, p = 0.04) (Fig. 3c, d). While molecular diversity declined (Fig. 3a, b) for arms and face, bacterial diversity (Fig. 3c, d) was less affected in the face and arms when participants did not use skin care products (T1–T3). The molecular diversity remained stable in the arms and face of female participants during common beauty products use (T4–T6) to immediately increase as soon as the volunteers went back to their normal routines (T7–T9) (WR test, p = 0.006 for the arms and face)(Fig. 3a, b). A higher molecular (Additional file 1: Figure S2A) and community (Additional file 1: Figure S2B) diversity was observed for armpits and feet of all individuals during the use of antiperspirant and foot powder (T4–T6) (WR test, molecular diversity: armpit p = 8.9e−33, foot p = 1.03e−11; bacterial diversity: armpit p = 2.14e−28, foot p = 1.26e−11), followed by a molecular and bacterial diversity decrease in the armpits when their regular personal beauty product use was resumed (T7–T9) (bacterial diversity: WR test, p = 4.780e−21, molecular diversity: WR test, p = 2.159e−21). Overall, our data show that refraining from using beauty products leads to lower molecular and bacterial diversity, while resuming the use increases their diversity. Distinct variations between male and female molecular and community richness were perceived at distinct body parts (Fig. 3a–d). Although the chemical diversity of personal beauty products does not explain these variations (Additional file 1: Figure S2C), differences observed between males and females may be attributed to many environmental and lifestyle factors including different original skin care and different frequency of use of beauty products (Additional file 2: Table S1), washing routines, and diet.

Longitudinal variation of skin metabolomics signatures

To gain insights into temporal metabolomics variation associated with beauty product use, chemical inventories collected over 9 weeks were subjected to multivariate analysis using the widely used Bray–Curtis dissimilarity metric (Fig. 4a–c, S3A). Throughout the 9-week period, distinct molecular signatures were associated to each specific body site: arm, armpit, face, and foot (Additional file 1: Figure S3A, Adonis test, p < 0.001, R² 0.12391). Mass spectrometric signatures displayed distinct individual trends at each specific body site (arm, armpit, face, and foot) over time, supported by their distinct locations in PCoA (principal coordinate analysis) space (Fig. 4a, b) and based on the Bray–Curtis distances between molecular profiles (Additional file 1: Figure S3B, WR test, all p values < 0.0001 from T0 through T9). This suggests a high molecular inter-individual variability over time despite similar changes in personal care routines. Significant differences in molecular patterns associated to ceasing (T1–T3) (Fig. 4b, Additional file 1: Figure S3C, WR test, T0 vs T1–T3 p < 0.001) and resuming the use of common beauty products (T4–T6) (Additional file 1: Figure S3C) were observed in the arm, face, and foot (Fig. 4b), although the armpit exhibited the most pronounced changes (Fig. 4b, Additional file 1: Figure S3D, E, random forest highlighting that 100% of samples from each phase were correctly predicted). Therefore, we focused our analysis on this region. Molecular changes were noticeable starting the first week (T1) of discontinuing beauty product use. As shown for armpits in Fig. 4c, these changes at the chemical level are specific to each individual, possibly due to the extremely personalized lifestyles before the study and match their original use of deodorant. Based on the initial use of underarm products (T0) (Additional file 2: Table S1), two groups of participants can be distinguished: a group of five volunteers who used stick deodorant as evidenced by the mass spectrometry data and another group of volunteers where we found few or no traces suggesting they never or infrequently used stick deodorants (Additional file 2: Table S1). Based on this criterion, the chemical trends shown in Fig. 4c highlight that individuals who used stick deodorant before the beginning of the study (volunteers 1, 2, 3, 9, and 12) displayed a more pronounced shift in their armpits’ chemistries as soon as they stopped using deodorant (T1–T3), compared to individuals who had low detectable levels of stick deodorant use (volunteers 4, 6, 7, and 10), or “rarely-to-never” (volunteers 5 and 11) use stick deodorants as confirmed by the volunteers (Additional file 1: Figure S3F, WR test, T0 vs T1–T3 all p values < 0.0001, with greater distance for the group of volunteers 1, 2, 3, 9, and 12, compared to volunteers 4, 5, 6, 7, 10, and 11). The most drastic shift in chemical profiles was observed during the transition period, when all participants applied the common antiperspirant on a daily basis (T4–T6) (Additional file 1: Figure S3D, E). Finally, the molecular profiles became gradually more similar to those collected before the experiment (T0) as soon as the participants resumed using their personal beauty products (T7–T9) (Additional file 1: Figure S3C), although traces of skin care products did last through the entire T7–T9 period in people who do not routinely apply these products (Fig. 4c).

Comparing chemistries detected in armpits at the end timepoints—when no products were used (T3) and during product use (T6)—revealed distinct molecular signatures characteristic of each phase (random forest highlighting that 100% of samples from each group were correctly predicted, see Additional file 1: Figure S3D, E). Because volunteers used the same antiperspirant during T4–T6, molecular profiles converged during that time despite individual patterns at T3 (Fig. 4b, c, Additional file 1: Figure S3D). These distinct chemical patterns reflect the significant impact of beauty products on skin molecular composition. Although these differences may in part be driven by beauty product ingredients detected on the skin (Additional file 1: Figure S1), we anticipated that additional host- and microbe-derived molecules may also be involved in these molecular changes.

To characterize the chemistries that vary over time, we used molecular networking, a MS visualization approach that evaluates the relationship between MS/MS spectra and compares them to reference MS/MS spectral libraries of known compounds [29, 30]. We recently showed that molecular networking can successfully organize large-scale mass spectrometry data collected from the human skin surface [18, 19]. Briefly, molecular networking uses the MScluster algorithm [31] to merge all identical spectra and then compares and aligns all unique pairs of MS/MS spectra based on their similarities where 1.0 indicates a perfect match. Similarities between MS/MS spectra are calculated using a similarity score, and are interpreted as molecular families [19, 24, 32–34]. Here, we used this method to compare and characterize chemistries found in armpits, arms, face, and foot of 11 participants. Based on MS/MS spectral similarities, chemistries highlighted through molecular networking (Additional file 1: Figure S4A) were associated with each body region with 8% of spectra found exclusively in the arms, 12% in the face, 14% in the armpits, and 2% in the foot, while 18% of the nodes were shared between all four body parts and the rest of spectra were shared between two body sites or more (Additional file 1: Figure S4B). Greater spectral similarities were highlighted between armpits, face, and arm (12%) followed by the arm and face (9%) (Additional file 1: Figure S4B).

Molecules were annotated with Global Natural Products Social Molecular Networking (GNPS) libraries [29], using accurate parent mass and MS/MS fragmentation patterns, according to level 2 or 3 of annotation defined by the 2007 metabolomics standards initiative [35]. Through annotations, molecular networking revealed that many compounds derived from steroids (Fig. 5a–d), bile acids (Additional file 1: Figure S5A-D), and acylcarnitines (Additional file 1: Figure S5E-F) were exclusively detected in the armpits. Using authentic standards, the identity of some pheromones and bile acids were validated to a level 1 identification with matched retention times (Additional file 1: Figure S6B, S7A, C, D). Other steroids and bile acids were either annotated using standards with identical MS/MS spectra but slightly different retention times (Additional file 1: Figure S6A) or annotated with MS/MS spectra match with reference MS/MS library spectra (Additional file 1: Figure S6C, D, S7B, S6E-G). These compounds were therefore classified as level 3 [35]. Acylcarnitines were annotated to a family of possible acylcarnitines (we therefore classify as level 3), as the positions of double bonds or cis vs trans configurations are unknown (Additional file 1: Figure S8A, B).

Among the steroid compounds, several molecular families were characterized: androsterone (Fig. 5a, b, d), androstadienedione (Fig. 5c), androstanedione (Additional file 1: Figure S6E), androstanolone (Additional file 1: Figure S6F), and androstenedione (Additional file 1: Figure S6G). While some steroids were detected in the armpits of several individuals, such as dehydroisoandrosterone sulfate (m/z 369.19, rt 247 s) (9 individuals) (Fig. 5a, Additional file 1: Figure S6A), androsterone sulfate (m/z 371.189, rt 261 s) (9 individuals) (Fig. 5b, Additional file 1: Figure S6C), and 5-alpha-androstane-3,17-dione (m/z 271.205, rt 249 s) (9 individuals) (Additional file 1: Figure S6E), other steroids including 1-dehydroandrostenedione (m/z 285.185, rt 273 s) (Fig. 5c, Additional file 1: Figure S6B), dehydroandrosterone (m/z 289.216, rt 303 s) (Fig. 5d, Additional file 1: Figure S6D), and 5-alpha-androstan-17.beta-ol-3-one (m/z 291.231, rt 318 s) (Additional file 1: Figure S6F) were only found in the armpits of volunteer 11 and 4-androstene-3,17-dione (m/z 287.200, rt 293 s) in the armpits of volunteer 11 and volunteer 5, both are male that never applied stick deodorants (Additional file 1: Figure S6G). Each molecular species exhibited a unique pattern over the 9-week period. The abundance of dehydroisoandrosterone sulfate (Fig. 5a, WR test, p < 0.01 for 7 individuals) and dehydroandrosterone (Fig. 5a, WR test, p = 0.00025) significantly increased during the use of antiperspirant (T4–T6), while androsterone sulfate (Fig. 5b) and 5-alpha-androstane-3,17-dione (Additional file 1: Figure S6E) display little variation over time. Unlike dehydroisoandrosterone sulfate (Fig. 5a) and dehydroandrosterone (Fig. 5d), steroids including 1-dehydroandrostenedione (Fig. 5c, WR test, p = 0.00024) and 4-androstene-3,17-dione (Additional file 1: Figure S6G, WR test, p = 0.00012) decreased in abundance during the 3 weeks of antiperspirant application (T4–T6) in armpits of male 11, and their abundance increased again when resuming the use of his normal skin care routines (T7–T9). Interestingly, even within the same individual 11, steroids were differently impacted by antiperspirant use as seen for 1-dehydroandrostenedione that decreased in abundance during T4–T6 (Fig. 5c, WR test, p = 0.00024), while dehydroandrosterone increased in abundance (Fig. 5d, WR test, p = 0.00025), and this increase was maintained during the last 3 weeks of the study (T7–T9).

In addition to steroids, many bile acids (Additional file 1: Figure S5A-D) and acylcarnitines (Additional file 1: Figure S5E-F) were detected on the skin of several individuals through the 9-week period. Unlike taurocholic acid found only on the face (Additional file 1: Figures S5A, S7A) and tauroursodeoxycholic acid detected in both armpits and arm samples (Additional file 1: Figures S5B, S7B), other primary bile acids such as glycocholic (Additional file 1: Figures S5C, S7C) and chenodeoxyglycocholic acid (Additional file 1: Figures S5D, S7D) were exclusively detected in the armpits. Similarly, acylcarnitines were also found either exclusively in the armpits (hexadecanoyl carnitines) (Additional file 1: Figures S5E, S8A) or in the armpits and face (tetradecenoyl carnitine) (Additional file 1: Figures S5F, S8B) and, just like the bile acids, they were also stably detected during the whole 9-week period.

Bacterial communities and their variation over time

Having demonstrated the impact of beauty products on the chemical makeup of the skin, we next tested the extent to which skin microbes are affected by personal care products. We assessed temporal variation of bacterial communities detected on the skin of healthy individuals by evaluating dissimilarities of bacterial collections over time using unweighted UniFrac distance [36] and community variation at each body site in association to beauty product use [3, 15, 37]. Unweighted metrics are used for beta diversity calculations because we are primarily concerned with changes in community membership rather than relative abundance. The reason for this is that skin microbiomes can fluctuate dramatically in relative abundance on shorter timescales than that assessed here. Longitudinal variations were revealed for the armpits (Fig. 6a) and feet microbiome by their overall trend in the PCoA plots (Fig. 6b), while the arm (Fig. 6c) and face (Fig. 6d) displayed relatively stable bacterial profiles over time. As shown in Fig. 6a–d, although the microbiome was site-specific, it varied more between individuals and this inter-individual variability was maintained over time despite same changes in personal care routine (WR test, all p values at all timepoints < 0.05, T5 p = 0.07), in agreement with previous findings that individual differences in the microbiome are large and stable over time [3, 4, 10, 37]. However, we show that shifts in the microbiome can be induced by changing hygiene routine and therefore skin chemistry. Changes associated with using beauty products (T4–T6) were more pronounced for the armpits (Fig. 6a, WR test, p = 1.61e−52) and feet (Fig. 6b, WR test, p = 6.15e−09), while little variations were observed for the face (Fig. 6d, WR test, p = 1.402.e−83) and none for the arms (Fig. 6c, WR test, p = 0.296).

A significant increase in abundance of Gram-negative bacteria including the phyla Proteobacteria and Bacteroidetes was noticeable for the armpits and feet of both females (Fig. 6e; Mann–Whitney U, p = 8.458e−07) and males (Fig. 6f; Mann–Whitney U, p = 0.0004) during the use of antiperspirant (T4–T6), while their abundance remained stable for the arms and face during that time (Fig. 6e, f; female arm p = 0.231; female face p value = 0.475; male arm p= 0.523;male face p = 6.848751e−07). These Gram-negative bacteria include Acinetobacter and Paracoccus genera that increased in abundance in both armpits and feet of females (Additional file 1: Figure S9A), while a decrease in abundance of Enhydrobacter was observed in the armpits of males (Additional file 1: Figure S9B). Cyanobacteria, potentially originating from plant material (Additional file 1: Figure S9C) also increased during beauty product use (T4–T6) especially in males, in the armpits and face of females (Fig. 6e) and males (Fig. 6f). Interestingly, although chloroplast sequences (which group phylogenetically within the cyanobacteria [38]) were only found in the facial cream (Additional file 1: Figure S9D), they were detected in other locations as well (Fig. 6e, f. S9E, F), highlighting that the application of a product in one region will likely affect other regions of the body. For example, when showering, a face lotion will drip down along the body and may be detected on the feet. Indeed, not only did the plant material from the cream reveal this but also the shampoo used for the study for which molecular signatures were readily detected on the feet as well (Additional file 1: Figure S10A). Minimal average changes were observed for Gram-positive organisms (Additional file 1: Figure S10B, C), although in some individuals the variation was greater than others (Additional file 1: Figure S10D, E) as discussed for specific Gram-positive taxa below.

At T0, the armpit’s microflora was dominated by Staphylococcus (26.24%, 25.11% of sequencing reads for females and 27.36% for males) and Corynebacterium genera (26.06%, 17.89% for females and 34.22% for males) (Fig. 7a—first plot from left and Additional file 1: Figure S10D, E). They are generally known as the dominant armpit microbiota and make up to 80% of the armpit microbiome [39, 40]. When no deodorants were used (T1–T3), an overall increase in relative abundance of Staphylococcus (37.71%, 46.78% for females and 30.47% for males) and Corynebacterium (31.88%, 16.50% for females and 44.15% for males) genera was noticeable (WR test, p < 3.071e−05) (Fig. 7a—first plot from left), while the genera Anaerococcus and Peptoniphilus decreased in relative abundance (WR test, p < 0.03644) (Fig. 7a—first plot from left and Additional file 1: Figure S10D, E). When volunteers started using antiperspirants (T4–T6), the relative abundance of Staphylococcus (37.71%, 46.78% females and 30.47% males, to 21.71%, 25.02% females and 19.25% males) and Corynebacterium (31.88%, 16.50% females and 44.15% males, to 15.83%, 10.76% females and 19.60% males) decreased (WR test, p < 3.071e−05) (Fig. 7a, Additional file 1: Figure S10D, E) and at the same time, the overall alpha diversity increased significantly (WR test, p = 3.47e−11) (Fig. 3c, d). The microbiota Anaerococcus (WR test, p = 0.0006018), Peptoniphilus (WR test, p = 0.008639), and Micrococcus (WR test, p = 0.0377) increased significantly in relative abundance, together with a lot of additional low-abundant species that lead to an increase in Shannon alpha diversity (Fig. 3c, d). When participants went back to normal personal care products (T7–T9), the underarm microbiome resembled the original underarm community of T0 (WR test, p = 0.7274) (Fig. 7a). Because armpit bacterial communities are person-specific (inter-individual variability: WR test, all p values at all timepoints < 0.05, besides T5 p n.s), variation in bacterial abundance upon antiperspirant use (T4–T6) differ between individuals and during the whole 9-week period (Fig. 7a—taxonomic plots per individual). For example, the underarm microbiome of male 5 exhibited a unique pattern, where Corynebacterium abundance decreased drastically during the use of antiperspirant (82.74 to 11.71%, WR test, p = 3.518e−05) while in the armpits of female 9 a huge decrease in Staphylococcus abundance was observed (Fig. 7a) (65.19 to 14.85%, WR test, p = 0.000113). Unlike other participants, during T0–T3, the armpits of individual 11 were uniquely characterized by the dominance of a sequence that matched most closely to the Enhydrobacter genera. The transition to antiperspirant use (T4–T6) induces the absence of Enhydrobacter (30.77 to 0.48%, WR test, p = 0.01528) along with an increase of Corynebacterium abundance (26.87 to 49.74%, WR test, p = 0.1123) (Fig. 7a—male 11).

In addition to the armpits, a decline in abundance of Staphylococcus and Corynebacterium was perceived during the use of the foot powder (46.93% and 17.36%, respectively) compared to when no beauty product was used (58.35% and 22.99%, respectively) (WR test, p = 9.653e−06 and p = 0.02032, respectively), while the abundance of low-abundant foot bacteria significantly increased such as Micrococcus (WR test, p = 1.552e−08), Anaerococcus (WR test, p = 3.522e−13), Streptococcus (WR test, p = 1.463e−06), Brevibacterium (WR test, p = 6.561e−05), Moraxellaceae (WR test, p = 0.0006719), and Acinetobacter (WR test, p = 0.001487), leading to a greater bacterial diversity compared to other phases of the study (Fig. 7b first plot from left, Additional file 1: Figure S10D, E, Fig. 3c, d).

We further evaluated the relationship between the two omics datasets by superimposing the principal coordinates calculated from metabolome and microbiome data (Procrustes analysis) (Additional file 1: Figure S11) [34, 41, 42]. Metabolomics data were more correlated with patterns observed in microbiome data in individual 3 (Additional file 1: Figure S11C, Mantel test, r = 0.23, p < 0.001), individual 5 (Additional file 1: Figure S11E, r = 0.42, p < 0.001), individual 9 (Additional file 1: Figure S11H, r = 0.24, p < 0.001), individual 10 (Additional file 1: Figure S11I, r = 0.38, p < 0.001), and individual 11 (Additional file 1: Figure S11J, r = 0.35, p < 0.001) when compared to other individuals 1, 2, 4, 6, 7, and 12 (Additional file 1: Figure S11A, B, D, F, G, K, respectively) (Mantel test, all r < 0.2, all p values < 0.002, for volunteer 2 p n.s). Furthermore, these correlations were individually affected by ceasing (T1–T3) or resuming the use of beauty products (T4–T6 and T7–T9) (Additional file 1: Figure S11A-K).

Overall, metabolomics–microbiome correlations were consistent over time for the arms, face, and feet although alterations were observed in the arms of volunteers 7 (Additional file 1: Figure S11G) and 10 (Additional file 1: Figure S11I) and the face of volunteer 7 (Additional file 1: Figure S11G) during product use (T4–T6). Molecular–bacterial correlations were mostly affected in the armpits during antiperspirant use (T4–T6), as seen for volunteers male 7 (Additional file 1: Figure S11G) and 11 (Additional file 1: Figure S11J) and females 2 (Additional file 1: Figure S11B), 9 (Additional file 1: Figure S11H), and 12 (Additional file 1: Figure S11K). This perturbation either persisted during the last 3 weeks (Additional file 1: Figure S11D, E, H, I, K) when individuals went back to their normal routine (T7–T9) or resembled the initial molecular–microbial correlation observed in T0 (Additional file 1: Figure S11C, G, J). These alterations in molecular–bacterial correlation are driven by metabolomics changes during antiperspirant use as revealed by metabolomics shifts on the PCoA space (Additional file 1: Figure S11), partially due to the deodorant’s chemicals (Additional file 1: Figure S1J, K) but also to changes observed in steroid levels in the armpits (Fig. 5A, C, D, Additional file 1: Figure S6G), suggesting metabolome-dependant changes of the skin microbiome. In agreement with previous findings that showed efficient biotransformation of steroids by Corynebacterium [43, 44], our correlation analysis associates specific steroids that were affected by antiperspirant use in the armpits of volunteer 11 (Fig. 5c, d, Additional file 1: Figure S6G) with microbes that may produce or process them: 1-dehydroandrostenedione, androstenedione, and dehydrosterone with Corynebacterium (r = − 0.674, p = 6e−05; r = 0.671, p = 7e−05; r = 0.834, p < 1e−05, respectively) (Additional file 1: Figure S12A, B, C, respectively) and Enhydrobacter (r = 0.683, p = 4e−05; r = 0.581, p = 0.00095; r = 0.755, p < 1e−05 respectively) (Additional file 1: Figure S12D, E, F, respectively).

Subject recruitment and sample collection

Twelve individuals between 25 and 40 years old were recruited to participate in this study, six females and six males. Female volunteer 8 dropped out of the study as she developed a skin irritation during the T1–T3 phase. All volunteers signed a written informed consent in accordance with the sampling procedure approved by the UCSD Institutional Review Board (Approval Number 161730). Volunteers were required to follow specific instructions during 9 weeks. They were asked to bring in samples of their personal care products they used prior to T0 so they could be sampled as well. Following the initial timepoint time 0 and during the first 3 weeks (week 1–week 3), volunteers were asked not to use any beauty products (Fig. 1b). During the next 3 weeks (week 4–week 6), four selected commercial beauty products provided to all volunteers were applied once a day at specific body part (deodorant for the armpits, soothing foot powder between the toes, sunscreen for the face, and moisturizer for front forearms) (Fig. 1b, Additional file 3: Table S2 Ingredient list of beauty products). During the first 6 weeks, volunteers were asked to shower with a head to toe shampoo. During the last 3 weeks (week 7–week 9), all volunteers went back to their normal routine and used the personal care products used before the beginning of the study (Fig. 1b). Volunteers were asked not to shower the day before sampling. Samples were collected by the same three researchers to ensure consistency in sampling and the area sampled. Researchers examined every subject together and collected metabolomics and microbiome samples from each location together. Samples were collected once a week (from day 0 to day 68—10 timepoints total) for volunteers 1, 2, 3, 4, 5, 6, 7, 9, 10, 11, and 12, and on day 0 and day 6 for volunteer 8. For individuals 4, 9, and 10, samples were collected twice a week. Samples collected for 11 volunteers during 10 timepoints: 11 volunteers × 10 timepoints × 4 samples × 4 body sites = 1760. Samples collected from 3 selected volunteers during 9 additional timepoints: 3 volunteers × 9 timepoints × 4 samples × 4 body sites = 432. All samples were collected following the same protocol described in [18]. Briefly, samples were collected over an area of 2 × 2 cm, using pre-moistened swabs in 50:50 ethanol/water solution for metabolomics analysis or in Tris-EDTA buffer for 16S rRNA sequencing. Four samples were collected from each body part right and left side. The locations sampled were the face—upper cheek bone and lower jaw, armpit—upper and lower area, arm—front of the elbow (antecubitis) and forearm (antebrachium), and feet—in between the first and second toe and third and fourth toe. Including personal care product references, a total of 2275 samples were collected over 9 weeks and were submitted to both metabolomics and microbial inventories.

Metabolite extraction and UPLC-Q-TOF mass spectrometry analysis

Skin swabs were extracted and analyzed using a previously validated workflow described in [18, 19]. All samples were extracted in 200 μl of 50:50 ethanol/water solution for 2 h on ice then overnight at − 20 °C. Swab sample extractions were dried down in a centrifugal evaporator then resuspended by vortexing and sonication in a 100 μl 50:50 ethanol/water solution containing two internal standards (fluconazole 1 μM and amitriptyline 1 μM). The ethanol/water extracts were then analyzed using a previously validated UPLC-MS/MS method [18, 19]. We used a ThermoScientific UltiMate 3000 UPLC system for liquid chromatography and a Maxis Q-TOF (Quadrupole-Time-of-Flight) mass spectrometer (Bruker Daltonics), controlled by the Otof Control and Hystar software packages (Bruker Daltonics) and equipped with ESI source. UPLC conditions of analysis are 1.7 μm C18 (50 × 2.1 mm) UHPLC Column (Phenomenex), column temperature 40 °C, flow rate 0.5 ml/min, mobile phase A 98% water/2% acetonitrile/0.1% formic acid (v/v), mobile phase B 98% acetonitrile/2% water/0.1% formic acid (v/v). A linear gradient was used for the chromatographic separation: 0–2 min 0–20% B, 2–8 min 20–99% B, 8–9 min 99–99% B, 9–10 min 0% B. Full-scan MS spectra (m/z 80–2000) were acquired in a data-dependant positive ion mode. Instrument parameters were set as follows: nebulizer gas (nitrogen) pressure 2 Bar, capillary voltage 4500 V, ion source temperature 180 °C, dry gas flow 9 l/min, and spectra rate acquisition 10 spectra/s. MS/MS fragmentation of 10 most intense selected ions per spectrum was performed using ramped collision induced dissociation energy, ranged from 10 to 50 eV to get diverse fragmentation patterns. MS/MS active exclusion was set after 4 spectra and released after 30 s.

Mass spectrometry data collected from the skin of 12 individuals can be found here MSV000081582.

LC-MS data processing

LC-MS raw data files were converted to mzXML format using Compass Data analysis software (Bruker Daltonics). MS1 features were selected for all LC-MS datasets collected from the skin of 12 individuals and blank samples (total 2275) using the open-source software MZmine [66]—see Additional file 4: Table S3 for parameters. Subsequent blank filtering, total ion current, and internal standard normalization were performed (Additional file 5: Table S4) for representation of relative abundance of molecular features (Fig. 2, Additional file 1: Figure S1), principal coordinate analysis (PCoA) (Fig. 4). For steroid compounds in Fig. 5a–d, bile acids (Additional file 1: Figure S5A-D), and acylcarnitines (Additional file 1: Figure S5E, F) compounds, crop filtering feature available in MZmine [66] was used to identify each feature separately in all LC-MS data collected from the skin of 12 individuals (see Additional file 4: Table S3 for crop filtering parameters and feature finding in Additional file 6: Table S5).

Heatmap in Fig. 2 was constructed from the bucket table generated from LC-MS1 features (Additional file 7: Table S6) and associated metadata (Additional file 8: Table S7) using the Calour command line available here: https://github.com/biocore/calour. Calour parameters were as follows: normalized read per sample 5000 and cluster feature minimum reads 50. Procrustes and Pearson correlation analyses in Additional file 1: Figures S10 and S11 were performed using the feature table in Additional file 9: Table S8, normalized using the probabilistic quotient normalization method [67].

16S rRNA amplicon sequencing

16S rRNA sequencing was performed following the Earth Microbiome Project protocols [68, 69], as described before [18]. Briefly, DNA was extracted using MoBio PowerMag Soil DNA Isolation Kit and the V4 region of the 16S rRNA gene was amplified using barcoded primers [70]. PCR was performed in triplicate for each sample, and V4 paired-end sequencing [70] was performed using Illumina HiSeq (La Jolla, CA). Raw sequence reads were demultiplexed and quality controlled using the defaults, as provided by QIIME 1.9.1 [71]. The primary OTU table was generated using Qiita (https://qiita.ucsd.edu/), using UCLUST (https://academic.oup.com/bioinformatics/article/26/19/2460/230188) closed-reference OTU picking method against GreenGenes 13.5 database [72]. Sequences can be found in EBI under accession number EBI: ERP104625 or in Qiita (qiita.ucsd.edu) under Study ID 10370. Resulting OTU tables were then rarefied to 10,000 sequences/sample for downstream analyses (Additional file 10 Table S9). See Additional file 11: Table S10 for read count per sample and Additional file 1: Figure S13 representing the samples that fall out with rarefaction at 10,000 threshold. The dataset includes 35 blank swab controls and 699 empty controls. The blank samples can be accessed through Qiita (qiita.ucsd.edu) as study ID 10370 and in EBI with accession number EBI: ERP104625. Blank samples can be found under the metadata category “sample_type” with the name “empty control” and “Swabblank.” These samples fell below the rarefaction threshold at 10,000 (Additional file 11: Table S10).

To rule out the possibility that personal care products themselves contained the microbes that induced the changes in the armpit and foot microbiomes that were observed in this study (Fig. 7), we subjected the common personal care products that were used in this study during T4–T6 also to 16S rRNA sequencing. The data revealed that within the limit of detectability of the current experiment, few 16S signatures were detected. One notable exception was the most dominant plant-originated bacteria chloroplast detected in the sunscreen lotion applied on the face (Additional file 1: Figure S9D), that was also detected on the face of individuals and at a lower level on their arms, sites where stable microbial communities were observed over time (Additional file 1: Figure S9E, F). This finding is in agreement with our previous data from the 3D cartographical skin maps that revealed the presence of co-localized chloroplast and lotion molecules [18]. Other low-abundant microbial signatures found in the sunscreen lotion include additional plant-associated bacteria: mitochondria [73], Bacillaceae [74, 75], Planococcaceae [76], and Ruminococcaceae family [77], but all these bacteria are not responsible for microbial changes associated to beauty product use, as they were poorly detected in the armpits and feet (Fig. 7).

To assess the origin of Cyanobacteria detected in skin samples, each Greengenes [72] 13_8 97% OTU table (per lane; obtained from Qiita [78] study 10,370) was filtered to only features with a p__Cyanobacteria phylum. The OTU maps for these tables—which relate each raw sequence to an OTU ID—were then filtered to only those observed p__Cyanobacteria OTU IDs. The filtered OTU map was used to extract the raw sequences into a single file. Separately, the unaligned Greengenes 13_8 99% representative sequences were filtered into two sets, first the set of representatives associated with c__Chloroplast (our interest database), and second the set of sequences associated with p__Cyanobacteria without the c__Chloroplast sequences (our background database). Platypus Conquistador [79] was then used to determine what reads were observed exclusively in the interest database and not in the background database. Of the 4,926,465 raw sequences associated with a p__Cyanobacteria classification (out of 318,686,615 total sequences), at the 95% sequence identity level with 100% alignment, 4,860,258 sequences exclusively recruit to full-length chloroplast 16S by BLAST [80] with the bulk recruiting to streptophytes (with Chlorophyta and Stramenopiles to a lesser extent). These sequences do not recruit non-chloroplast Cyanobacteria full length 16S.

Half-life calculation for metabolomics data

In order to estimate the biological half-life of molecules detected in the skin, the first four timepoints of the study (T0, T1, T2, T3) were considered for the calculation to allow the monitoring of personal beauty products used at T0. The IUPAC’s definition of biological half-life as the time required to a substance in a biological system to be reduced to half of its value, assuming an approximately exponential removal [81] was used. The exponential removal can be described as C_(t) = C₀e^−tλ where t represents the time in weeks, C₀ represents the initial concentration of the molecule, C_(t) represents the concentration of the molecule at time t, and λ is the rate of removal [http://onlinelibrary.wiley.com/doi/10.1002/9780470140451.ch2/summary]. The parameter λ was estimated by a mixed linear effects model in order to account for the paired sample structure. The regression model tests the null hypothesis that λ is equal to zero and only the significant (p value < 0.05) parameters were considered.

Principal coordinate analysis

We performed principal coordinate analysis (PCoA) on both metabolomics and microbiome data. For metabolomics, we used MS1 features (Additional file 5: Table S4) and calculated Bray–Curtis dissimilarity metric using ClusterApp (https://github.com/mwang87/q2_metabolomics).

For microbiome data, we used rarefied OTU table (Additional file 10: Table S9) and used unweighted UniFrac metric [36] to calculate beta diversity distance matrix using QIIME2 (https://qiime2.org). Results from both data sources were visualized using Emperor (https://biocore.github.io/emperor/) [28].

Molecular networking

Molecular networking was generated from LC-MS/MS data collected from skin samples of 11 individuals MSV000081582, using the Global Natural Products Social Molecular Networking platform (GNPS) [29]. Molecular network parameters for MS/MS data collected from all body parts of 11 individuals during T0–T9 MSV000081582 are accessible here http://gnps.ucsd.edu/ProteoSAFe/status.jsp?task=284fc383e4c44c4db48912f01905f9c5. Molecular network parameters for MS/MS data collected from armpits T0–T3 MSV000081582 and deodorant used by individual 1 and 3 MSV000081580 can be found here http://gnps.ucsd.edu/ProteoSAFe/status.jsp?task=f5325c3b278a46b29e8860ec57915ad and here http://gnps.ucsd.edu/ProteoSAFe/status.jsp?task=aaa1af68099d4c1a87e9a09f398fe253, respectively. Molecular networks were exported and visualized in Cytoscape 3.4.0. [82]. Molecular networking parameters were set as follows: parent mass tolerance 1 Da, MS/MS fragment ion tolerance 0.5 Da, and cosine threshold 0.65 or greater, and only MS/MS spectral pairs with at least 4 matched fragment ions were included. Each MS/MS spectrum was only allowed to connect to its top 10 scoring matches, resulting in a maximum of 10 connections per node. The maximum size of connected components allowed in the network was 600, and the minimum number of spectra required in a cluster was 3. Venn diagrams were generated from Cytoscape data http://gnps.ucsd.edu/ProteoSAFe/status.jsp?task=284fc383e4c44c4db48912f01905f9c5 using Cytoscape [82] Venn diagram app available here http://apps.cytoscape.org/apps/all.

Shannon molecular and bacterial diversity

The diversity analysis was performed separately for 16S rRNA data and LC-MS data. For each sample in each feature table (LC-MS data and microbiome data), we calculated the value of the Shannon diversity index. For LC-MS data, we used the full MZmine feature table (Additional file 5: Table S4). For microbiome data, we used the closed-reference BIOM table rarefied to 10,000 sequences/sample. For diversity changes between timepoints, we aggregated Shannon diversity values across groups of individuals (all, females, males) and calculated mean values and standard errors. All successfully processed samples (detected features in LC-MS or successful sequencing with 10,000 or more sequences/sample) were considered.

Beauty products and chemical standards

Samples (10 mg) from personal care products used during T0 and T7–T9 MSV000081580 (Additional file 2: Table S1) and common beauty products used during T4–T6 MSV000081581 (Additional file 3: Table S2) were extracted in 1 ml 50:50 ethanol/water. Sample extractions were subjected to the same UPLC-Q-TOF MS method used to analyze skin samples and described above in the section “Metabolite extraction and UPLC-Q-TOF mass spectrometry analysis.” Authentic chemical standards MSV000081583 including 1-dehydroandrostenedion (5 μM), chenodeoxyglycocholic acid (5 μM), dehydroisoandrosterone sulfate (100 μM), glycocholic acid (5 μM), and taurocholic acid (5 μM) were analyzed using the same mass spectrometry workflow used to run skin and beauty product samples.

Monitoring beauty product ingredients in skin samples

In order to monitor beauty product ingredients used during T4–T6, we selected only molecular features present in each beauty product sample (antiperspirant, facial lotion, body moisturizer, soothing powder) and then filtered the aligned MZmine feature table (Additional file 5: Table S4) for the specific feature in specific body part samples. After feature filtering, we selected all features that had a higher average intensity on beauty product phase (T4–T6) compared to non-beauty product phase (T1–T3). The selected features were annotated using GNPS dereplication output http://gnps.ucsd.edu/ProteoSAFe/status.jsp?task=69319caf219642a5a6748a3aba8914df, plotted using R package ggplot2 (https://cran.r-project.org/web/packages/ggplot2/index.html) and visually inspected for meaningful patterns.

Random forest analysis

Random forest analysis was performed in MetaboAnalyst 3.0 online platform http://www.metaboanalyst.ca/faces/home.xhtml. Using LC-MS1 features found in armpit samples collected on T3 and T6. Random forest parameters were set as follows: top 1000 most abundant features, number of predictors to try for each node 7, estimate of error rate (0.0%).

BugBase analysis

To determine the functional potential of microbial communities within our samples, we used BugBase [83]. Because we do not have direct access to all of the gene information due to the use of 16S rRNA marker gene sequencing, we can only rely on phylogenetic information inferred from OTUs. BugBase takes advantage of this information to predict microbial phenotypes by associating OTUs with gene content using PICRUSt [84]. Thus, using BugBase, we can predict such phenotypes as Gram staining, or oxidative stress tolerance at each timepoint or each phase. All statistical analyses in BugBase are performed using non-parametric differentiation tests (Mann–Whitney U).

Taxonomic plots

Rarefied OTU counts were collapsed according to the OTU’s assigned family and genus name per sample, with a single exception for the class of chloroplasts. Relative abundances of each family-genus group are obtained by dividing by overall reads per sample, i.e., 10,000. Samples are grouped by volunteer, body site, and time/phase. Abundances are aggregated by taking the mean overall samples, and resulting abundances are again normalized to add up to 1. Low-abundant taxa are not listed in the legend and plotted in grayscale. Open-source code is available at https://github.com/sjanssen2/ggmap/blob/master/ggmap/snippets.py

Dissimilarity-based analysis

Pairwise dissimilarity matrices were generated for metabolomics and 16S metagenomics quantification tables, described above, using Bray–Curtis dissimilarity through QIIME 1.9.1 [71]. Those distance matrices were used to perform Procrustes analysis (QIIME 1.9.1), and Mantel test (scikit-bio version 0.5.1) to measure the correlation between the metabolome and microbiome over time. The metabolomics dissimilarities were used to perform the PERMANOVA test to assess the significance of body part grouping. The PCoA and Procrustes plots were visualized in EMPeror. The dissimilarity matrices were also used to perform distance tests, comparing the distances within and between individuals and distances from time 0 to times 1, 2, and 3 using Wilcoxon rank-sum tests (SciPy version 0.19.1) [19].

Statistical analysis for molecular and microbial data

Statistical analyses were performed in R and Python (R Core Team 2018). Monotonic relationships between two variables were tested using non-parametric Spearman correlation tests. The p values for correlation significance were subsequently corrected using Benjamini and Hochberg false discovery rate control method. The relationship between two groups was tested using non-parametric Wilcoxon rank-sum tests. The relationship between multiple groups was tested using non-parametric Kruskal–Wallis test. The significance level was set to 5%, unless otherwise mentioned, and all tests were performed as two-sided tests.