Skin metatranscriptomics reveals landscape of variation in microbial activity and gene expression across the human body

18 The skin microbiome plays an important role in immune homeostasis and skin health, 19 and yet our understanding of in vivo microbial gene activity is hindered by the lack of 20 a robust, non -invasive protocol for metatranscriptomics across skin sites . 21 Circumventing the challenges of low microbial biomass, host contamination, and RNA 22 stability, we developed a clinically tractable skin metatranscriptomics workflow that 23 provides high technical reproducibility of profiles (Pearson r>0.95), uniform coverage 24 across gene bodies, and strong enrichment of microbial mRNAs ( 2.5-40). Applying 25 this protocol to a cohort of healthy adults (n=27) across five different skin sites (n=102, 26 paired metatranscriptomes and metagenomes) , identified a striking di vergence 27 between transcriptomic and genomic abundances, with Staphylococcus species and 28 the skin fungi Malassezia having an outsized contribution to the metatranscriptomic 29 landscape at most sites despite their modest representation in metagenomes. 30 Species-level analysis showed skin site-specific enrichment of gene expression (e.g. 31 increased levels of secreted fungal phospholipase C on cheeks relative to scalp), and 32 revealed how key pathways were transcriptionally active in vivo (e.g. propionate and 33 4-aminobutyrate metabolism, potentially impacting skin barrier function ). Gene-level 34 analysis identified diverse antimicrobial genes transcribed by skin commensals in situ, 35 including several uncharacterized bacteriocins, some of which are expressed at levels 36 comparable to known antimicrobial genes . Correlation of microbial gene expression 37 with organismal abundances uncovered >20 gene s that putatively mediate 38 interactions between microbes (e.g. a secreted Malassezia restricta protein with 39 strongly negative in vivo association with Cutibacterium acnes; Spearman ρ>0.7). This 40 work showcases the potential for leveraging skin metatranscriptomics to identify 41 microbes whose activities play an outsized role in the community, and for uncovering 42 pivotal microbial pathways and biomarkers linked to skin health and disease. 43 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint

44 The human skin is home to diverse communities of microorganisms (bacteria, fungi, 45 and viruses) that can interact with each other and the host to impact the skin 46 microenvironment, immune homeostasis and skin health 1. In recent year s, o ur 47 understanding of the role of the skin microbiome in various diseases has greatly 48 benefited from the increasing accessibility of whole metagenome sequencing to 49 identify key organisms and associated genetic potential differences2–4. This has been 50 complemented by a growing number of ex vivo5,6 and in vitro7,8 studies that provide 51 insights into mechanistic pathways through which the crosstalk between the host and 52 the microbiome may be mediated . However, our understanding of whether these 53 pathways are indeed utilized in vivo remains limited, as metagenomics only estimates 54 functional potential of microbes 9 (DNA content) . Metagenomic DNA signals are a 55 composite from living and dead cells , and even among living microbes, their genes 56 can be variably expressed or transcriptionally silent in response to nutrient availability 57 and other environmental cues 10–12. Consequently, much remains unknown about 58 which skin microbial species and pathways are transcriptionally active in vivo, their 59 variability across skin sites and subjects, and their relationship to the skin 60 metagenome. 61 Metatranscriptomics, which assays the pool of messenger RNAs (mRNAs) in a 62 microbial community, has been used to study the transcriptional activity of microbiota 63 in environments as diverse as the gut 13 and ocean water 14, but has not yet been 64 broadly applied to the skin microbiome. Skin metatranscriptomics is hampered by the 65 lack of a robust, non -invasive protocol that can accommodate a range of skin sites 66 which have low microbial biomass , but substantial host and environmental 67 contamination. This is because human skin is relatively sparsely colonized by 68 microbes, with an estimated average density of 10 3-104 prokaryotes/cm2, several 69 orders of magnitude lower than that of the gut 15. Till date, there has been only one 70 other published study which has leveraged RNA sequencing to reveal how microbes 71 contribute to acne, albeit restricted to analyzing data from a single microbe 72 (Cutibacterium acnes ), and specifically applied for nose skin follicles using pore -73 stripping16. While biopsies are not limited to follicles and can capture more microbial 74 biomass17, their invasive nature makes them impractical for adoption in large -scale 75 clinical studies . Consequently, data from transcriptionally active microbes and 76 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint microbial pathways across different skin sites is scarce and highlights the need for a 77 generalized protocol to characterize skin metatranscriptomes18. 78 To address these challenges, we developed a robust workflow for non-invasive 79 sampling and skin metatranscriptomics across body sites that demonstrates high 80 technical reproducibility (Pearson r>0.95), uniform coverage across bacterial and 81 fungal gene bodies, and strong enrichment of microbial mRNAs (2.5-40). We further 82 developed a data analysis workflow with rigorous control of “kitome” contaminants and 83 taxonomic misclassification artifacts. The workflow is customized for profiling the skin 84 metatranscriptome with the use of an updated skin microbial gene catalog based on 85 metagenome-assembled genomes (MAGs) for diverse populations19. Leveraging this 86 capability, we present the first multi-site metatranscriptomic survey of healthy human 87 skin (n=27 subjects) from physiologically diverse skin sites (n=5; scalp, cheek, volar 88 forearm, antecubital fossae and toe web ). Analysis of paired metagenomes and 89 metatranscriptomes (n=102 paired libraries; n=260 libraries in total) revealed a striking 90 divergence between transcriptomic and genomic abundances, with Staphylococcus 91 species and the skin fungi Malassezia having an outsized contribution to the 92 metatranscriptomic landscape at most sites despite their l imited representation in 93 metagenomes. Species-level analysis showed skin site -specific enrichment of gene 94 expression (e.g. increased levels of secreted fungal phospholipase C on cheeks 95 relative to scalp), as well as revealed transcriptional activity in vivo of key pathways 96 (e.g. propionate and 4 -aminobutyrate metabolism, potentially impacting skin barrier 97 function). Gene-level analysis identified diverse antimicrobial genes transcribed by 98 skin commensals in situ , including several uncharacterized bacteriocins, some of 99 which are expressed at levels comparable to known microbially produced antimicrobial 100 genes. Correlation of microbial gene expression with organismal abundances 101 uncovered >20 genes that putatively mediate interactions between microbes (e.g. a 102 secreted Malassezia restricta protein with strongly negative in vivo association with 103 Cutibacterium acnes; Spearman ρ>0.7). Overall, this work highlights the importance 104 of metatranscriptomics for a holistic view of active species, expressed microbial 105 functions, key pathways and microbe-microbe interactions occurring in situ on skin. 106 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint

107 Development of a robust skin metatranscriptomics workflow 108 We optimized an experimental workflow that is robust for metatranscriptomics across 109 different skin sites by systematically and extensively testing different sampling tools 110 (skin tapes versus swabs), lysis conditions and RNA purification techniques 111 (Supplementary Data 1; Methods). Our final optimized protocol utilized skin swabs, 112 sample collection in DNA/RNA shield, bead beating, rRNA depletion using a custom 113 oligonucleotide mix and a direct-to-column TRIzol purification step ( Methods). The 114 robustness of our protocol was assessed with a pilot cohort consisting of both a cross-115 sectional group (n=2 4 samples) and a longitudinal group sampled across three 116 consecutive days (n=45 samples), representing five distinct skin microenvironments: 117 scalp (Sc), cheek (Ch), volar forearm (Vf), antecubital fossae (Ac) and toe web (Tw). 118 A notably high proportion of metatranscriptomic libraries were successfully sequenced 119 in this cohort (66/69=95%), enabling generation of a target million microbial reads per 120 sample in most samples (>84%, median=2.2 million microbial reads, 0.66 Gbp; 121 Supplementary Data 1; Methods). Depletion of rRNAs resulted in substantial 122 enrichment (2.5-40) of reads corresponding to non-rRNAs (median>79.5% of reads; 123 Figure 1A ). Analysis of the resulting reads showed that the corresponding s kin 124 metatranscriptomes were highly reproducible across different skin site s, with very 125 strong consistency of species profiles ( Sorensen similarity≥0.98; Figure 1B) as well 126 as microbial gene expression (Pearson’s r≥0.99; Figure 1 C) between technical 127 replicates. Analysis of samples from the longitudinal group highlighted that 128 metatranscriptomes can exhibit substantial temporal stability at the gene level (median 129 Pearson’s r≥0.897 within individuals; Figure 1C), while being slightly more variable at 130 the species level (median Sorensen similarity≥ 0.768; Figure 1B), with the temporal 131 gene/species level variation being significantly lower than inter-individual variation at 132 the same skin site (Wilcoxon p-value<10-11; Figure 1B-C). 133 We next developed a computational workflow that could annotate skin 134 metatranscriptomic reads with high sensitivity, while accounting for potential 135 contamination signals and off -target matches. We firstly noted that using a skin-136 specific microbial gene catalog19 (IHSMGC) and a custom workflow, as opposed to a 137 widely used general-purpose workflow (HUMAnN320), resulted in a significantly higher 138 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint median percentage of reads being functionally annotated (81% vs 60%, Wilcoxon p-139 value<3.110-5; Supplementary Figure 1A ; Methods). We then used data from 140 negative controls as well as prior reports 21 to systematically identify potential DNA 141 contaminant taxa22 (i.e. the kitome) and filtered corresponding reads from our analysis 142 (Supplementary Figure 2A; Methods). Correlation analysis confirmed that skin 143 microbes and potential kitome taxa (e.g. Achromobacter, Bradyrhizobium, 144 Mycolibacterium, Mycobacterium and Brevundimonas species) formed distinct 145 clusters of co-occurring taxa (Supplementary Figure 2B). To account for potential 146 taxonomic classification errors, particularly in low complexity or misassembled regions 147 of microbial genomes 23, we compared a measure of unique matches in the genome 148 based on minimizers with the total read count 23,24 to identify false-positive taxa 149 (Methods). Using data from spike-ins we found that an empirically determined 150 threshold of unique minimizers per million microbial reads could discriminate false 151 from true positive taxa at relative abundances as low as 0.1%, over a range of read 152 counts (10 4-106 reads; Supplementary Figure 2C; Methods) and was therefore 153 consistently applied as a filter. 154 We applied the combined experimental and computational skin 155 metatranscriptomics workflow to the full cohort of 27 healthy individuals at 5 skin sites 156 collected specifically for this study to provide comprehensive in vivo characterization 157 of microbial gene expression on skin (n=135 paired samples for metagenomics and 158 metatranscriptomics; Figure 1D, Supplementary Data 1; Methods). In the full cohort, 159 the success rate for metatranscriptomic ( 102/135=75%) and metagenomic 160 (130/135=96%) libraries was moderate-to-high, emphasizing the robustness of this 161 protocol across a group of individuals. Typically, >1Gbp of de-duplicated non-rRNA 162 sequencing data (median of 3.7 million read pairs) was generated per library similar 163 to other metatranscriptomic studies 13,25, with a relatively high median RNA quality 164 across skin sit es (DV200≥76; Figure 1D, Supplementary Figure 1B ). In addition, 165 paired metagenomes were sequenced to sufficient depth to obtain a median of 7 166 million microbial read pairs after filtering for human reads (Supplementary Figure 1C). 167 Rarefaction analysis confirmed that libraries in the full cohort were typically sequenced 168 at sufficient depths for representing active microbial functions/orthologous groups (>1 169 million read pairs ; Supplementary Figure 3). Interestingly, t he proportion of non -170 human reads was found to be significantly higher in metatranscriptomes versus 171 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint metagenomes (98% vs 10%, Wilcoxon signed rank p -value<0.05; Supplementary 172 Figure 1 D), underscoring the feasibility of skin metatranscriptomic sequencing . 173 Microbial rRNAs across all the sites were effectively depleted (2 -25% of remaining 174 RNA compared to 80 -90% in a typical cell 26) during the library preparation process 175 (Figure 1E). In addition, the sequenced reads exhibited relatively even coverage 176 across bacterial (Figure 1F) as well as fungal (Figure 1G) gene bodies in all five skin 177 sites, and across a range of metatranscriptomic read depths. Similar to the pilot cohort, 178 annotation rates for reads were found to be moderate -to-high across various body 179 sites (median=69-80%, Supplementary Figure 1E). In addition, most of the reads 180 mapping to the genome of the common skin fungi Malassezia globosa were in coding 181 regions of mRNAs (>80%), with a minority mapping to intergenic regions and introns, 182 indicating that our metatranscriptome libraries achieved DNA depletion levels similar 183 to other high-quality M. globosa RNA-seq datasets derived from in vitro cultures27 184 (Figure 1 H). Overall, these results emphasize that the generalized skin 185 metatranscriptomics protocol presented here enables robust and reproducible profiling 186 of microbial mRNAs across a wide spectrum of body sites and subjects. 187 Skin metatranscriptomics identifies niche-specific active species and functions 188 distinct from metagenomes 189 Different skin microenvironments (e.g. sebaceous, moist, or dry) greatly influence the 190 composition of microbes that are present4, but much remains unknown about species 191 or gene activities in vivo. While some studies have reported discordance between RNA 192 and DNA abundances for specific species and gene families in gut and environmental 193 microbiomes14,25, the scale and extent to which gene abundances and transcript levels 194 correlate in the skin microbiome is not known. In this context, we observed a striking 195 disparity between the most active species in skin metatranscriptomes versus the most 196 highly abundant species in skin metagenomes (Figure 2A). For example, 197 Cutibacterium acnes which is a dominant component of the metagenome at most skin 198 sites (46-90% median relative abundance; except toe webs) has a relatively modest 199 contribution to the metatranscriptome (2-31% median relative abundance). In contrast, 200 while the skin fungi Malassezia restricta and Malassezia globosa were present at 201 relatively low metagenomic abundances relative to bacteria (3-8% and 0. 1-12% 202 median relative abundance , respectively), they contribute substantially to 203 metatranscriptomes in a niche dependent manner . M. restricta RNAs are heavily 204 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint represented in cheek and scalp sites (23-30% median relative abundance), whereas 205 M. globosa RNAs predominate on the scalp, antecubital fossae and volar forearms 206 (21-81% median relative abundance). Toe webs were the exception to this observation, 207 with a distinct microbial composition dominated by Staphylococcus hominis and 208 Staphylococcus epidermidis in metagenomes and metatranscriptomes ( Figure 2A). 209 While n ormalized RNA counts of many commonly reported skin microbes varied 210 positively with genomic DNA counts , several species had relatively large differences 211 (≥4-fold) in RNA and DNA abundances, with disproportionately higher contributions to 212 transcriptional activity (Supplementary Figure 4). M. restricta and M. globosa’s 213 outsized contribution to the active biomass across different skin sites is likely driven 214 by the larger cell volume of eukaryotes compared to prokaryotes 28. Interestingly, 215 Staphylococcus and Corynebacterium species also contributed proportionally more 216 RNAs at the scalp, cheek, volar forearm and antecubital fossa a fter accounting for 217 their metagenomic abundances (Supplementary Figure 4). The opposite was true 218 for C. acnes and Micrococcus luteus, likely reflecting a relatively low proportion of cells 219 present that are transcriptionally active . Overall, the stark differences between DNA 220 and RNA abundances for various skin microbes underlines the need for skin 221 metatranscriptomic measurements to characterize in vivo microbial activity. 222 To profile the landscape of skin microbial activity across different niches, we 223 categorized skin microbial species in terms of their prevalence as belonging to ‘core’ 224 (present in >75% of samples), ‘common’ (present in 50-75% of samples) and ‘variable’ 225 (present in <50% of samples) components of metagenomes and metatranscriptomes 226 at each site. The median ratio of RNA/DNA levels per species was used to quantify 227 transcriptional activity while accounting for variations in genom ic abundances and 228 sequencing depth ( Figure 2B ). Several species were found to be core in the 229 metagenomes of all or most surveyed skin sites such as Corynebacterium 230 kefirresidentii, C. acnes, M. globosa, M. restricta, Staphyloccocus capitis and S. 231 epidermidis (Figure 2B). However fewer species were core components of the skin 232 metatranscriptome across sites (e.g. C. acnes , M. globosa and M. restricta ), with 233 Staphylococcus species in particular exhibiting site -specific activity (e.g. S. 234 epidermidis on cheek with median activity >17, S. hominis on volar forearm with 235 median activity >9.5, and S. capitis on scalp with median activity > 3.5), and variable 236 prevalence across sites (e.g. in antecubital fossae; Figure 2B). 237 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint The same species also exhibited site-specific differences in overall gene 238 expression activity that could reflect their preferences for certain skin 239 microenvironments. Microbes that were prevalently found across multiple skin sites 240 such as C. acnes, S. capitis and S. epidermidis showed higher overall transcriptional 241 activity in sebaceous ( scalp and cheek) versus non -sebaceous (antecubital fossae, 242 volar forearm and toe web) sites, reflecting a general preference for environments 243 richer in host lipids 29,30 (Supplementary Figure 5). Intriguingly, although M. globosa 244 and M. restricta are closely related lipophilic fungi, the latter showed larger variation in 245 transcriptional activity between sebaceous and non -sebaceous sites, indicating that 246 M. restricta may have greater sensitivity to host lipid availability (Figure 2B , 247 Supplementary Figure 6). This is also in line with other reports demonstrating that M. 248 globosa can be readily isolated from facial skin, the upper trunk and arms, while M. 249 restricta recovery is mostly restricted to the scalp and face 31. Within skin sites, 250 individual heterogeneity is also evident in the species that are dominant in the 251 metatranscriptome. For instance, while transcripts from Malassezia spp. make up 252 most of the metatranscriptomes of the antecubital fossae (median relative abundance 253 98%), two individuals exhibited metatranscriptomes primarily dominated by S. hominis 254 and C. acnes respectively (SMT003_Ac and SMT011_Ac; Supplementary Figure 7, 255 3rd and 8th column from the left). These observations highlight that the distribution of 256 active skin microbes is shaped by skin site-specific differences but also individual-257 specific variations in the microenvironment. 258 To determine what functions are relatively important for microbial communities 259 in distinct niches, we performed differential expression analysis of microbial genes 260 clustered into orthologous groups (OGs), accounting for variations in DNA 261 abundances (Methods). Compared to those on a sebum rich environment such as the 262 cheek, bacteria colonizing the relatively dry volar forearms were found to upregulate 263 gene clusters important for glucose catabolism, energy generation and protein 264 synthesis, reflecting differences in resource availability (Supplementary Figure 8A, 265 Supplementary Figure 9A). Malassezia fungi colonizing the cheeks upregulated 266 gene clusters involved in mitotic growth, aromatic compound biosynthesis and protein 267 modification relative to those in the volar forearms , consistent with increased lipid 268 availability on the cheek for fungal growth and metabolism 28,32 (Supplementary 269 Figure 8A, Supplementary Figure 9A). Gene clusters involved in the citric acid cycle, 270 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint amino acid metabolism, heme biosynthesis, RNA modification and kinase activities 271 were upregulated in the toe webs relative to the volar forearms, indicative of microbes 272 adapting their metabolism to an amino-acid rich environment there (Supplementary 273 Figure 8B). Amino acid metabolic pathways are noteworthy because s weat is a rich 274 source of free amino acids33 and the toe webs represent an environment high in sweat 275 content relative to the volar forearm. The biosynthesis of heme, an important co-factor 276 for many microbial processes involved in Staphylococcal colonization is also 277 associated with increased availability of amino acids 34,35. Consequently, most 278 enzymes in the heme biosynthetic pathway from glutamate were upregulated by two-279 fold or more in toe webs compared to the volar forearm ( Supplementary Figure 9B-280 C). Consistent with the more exposed nature of volar forearms and lower 281 concentrations of protective secretions compared to cheeks 36, microbes on volar 282 forearms had upregulated expression of genes important for antioxidant protection 283 relative to cheeks (Supplementary Figure 8A) and toe webs (Supplementary Figure 284 8B). These examples highlight the utility of metatranscriptomics for identifying key 285 niche-specific functions in the skin microbiome. 286 We evaluated the importance of different species to core (present in >75% of 287 individuals) biochemical pathways across different skin sites, as well as within-sample 288 (alpha) and between-sample (beta) contributional diversities per pathway9 (Methods). 289 At all sites besides the toe webs, Malassezia fungi were the predominant (>50% 290 contribution) effectors of a substantial fraction (median 56%) of all core pathways 291 compared to other bacterial genera (Wilcoxon p -value<0.05; Supplementary Figure 292 10A). Within bacteria genera, Staphylococcus was the major contributor to more core 293 pathways than Cutibacterium at moist environments such as antecubital fossae and 294 toe webs (Wilcoxon p -value<0.05; Supplementary Figure 10A ). At any given site, 295 some pathways might be redundantly expressed by multiple microbes while other s 296 might only be driven by one or a few species. Most core microbial pathways on skin 297 were expressed by a few species within the same individual, but the same functions 298 were expressed by diverse sets of species across individuals, indicating high 299 functional plasticity in skin microbial communities (Figure 2C). For example, 300 staphyloferrin A biosynthesis from L -ornithine is usually expressed by 1-2 301 Staphylococcus species in individual toe web communities but multiple possible 302 Staphylococcus species can drive expression of this pathway across individuals 303 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint (Figure 2D ). Some metabolic functions were solely represented by fungi such as 304 pathways for beta oxidation of very long chain fatty acids (VLCFA ) being primarily 305 expressed by Malassezia species on both sebaceous (scalp) and non-sebaceous sites 306 (antecubital fossae) within and across individuals (Supplementary Figure 10B). 307 Interestingly, other metabolic functions such as galactose degradation or arginine 308 biosynthesis could be shared between skin bacteria and fungi, especially in nutrient -309 rich niches like the scalp or cheeks ( Supplementary Figure 10C-D). While bacterial 310 arginine has been described as a source of natural moisturizing factor (NMFs) for skin 311 health37, our dataset indicates that fungi may be substantial contributors to the arginine 312 pool. Overall, these results showcase how skin metatranscriptomics can be a 313 complementary approach to metagenomics and demonstrates how different genera or 314 species can act as the predominant effectors of niche -specific metabolic functions 315 across skin sites. 316 Species-level transcriptome analysis identifies signatures of microbial adaption 317 to in vivo nutrient availability 318 Given the presence of niche-specific metatranscriptomic signatures, we next 319 conducted species -level differential gene expression analyses to identify 320 transcriptional variation in different environments that may highlight specific metabolic 321 pathways essential for supporting in vivo colonization. We focused on RNA sequences 322 belonging to the skin commensals M. restricta, M. globosa, C . acnes and S. 323 epidermidis, because of their abundant transcript levels and importance for skin 324 health28,38. In particular, Malassezia are highly represented in the metatranscriptomes 325 of sebaceous sites (scalp and cheeks), while S. epidermidis can colonize both 326 sebaceous and non-sebaceous sites (especially toe webs) and can be readily cultured 327 in vitro . Differential expression analysis of M. restricta transcripts on scalp versus 328 cheek showed that fungal genes and pathways associated with fructose and mannose 329 metabolism, metabolism in diverse environments (e.g. ACU7_1 and acuE of the 330 glyoxylate cycle; FBA1 and KGD1 of the citric acid cycle; TLK1 of the pentose 331 phosphate pathway), and free fatty acid breakdown in peroxisomes (e.g. POX2, POT1, 332 PEX11B and ACOT8) were upregulated on scalp relative to cheeks ( Figure 3A-B, 333 Supplementary Data 2). Conversely, genes and pathways associated with the 334 breakdown of glycerides such as ether lipid and glycerophospholipid metabolism, as 335 well as various secreted phospholipase C enzymes were relatively upregulated on 336 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint cheeks (PLC_1-7; Figure 3A-B, Supplementary Data 2). This is consistent with 337 distinct metaboli te and lipid profiles on scalp and cheeks , with the former being 338 relatively richer in nutritionally complex apocrine secretions and free fatty acids, while 339 the latter being richer in glycerides 39,40. Intriguingly, although M. globosa is another 340 lipophilic skin fungi with comparable metatranscriptomic abundances to M. restricta on 341 sebaceous sites (Figure 2A), there were few differences in gene expression between 342 scalp and cheek sites for M. globosa (<30; Figure 3C). Another core skin commensal, 343 C. acnes, had a similar number of differentially expressed features between the scalp 344 and cheeks as M. restricta, but showed fewer significantly enriched pathways, with 345 only those related to ribosome biogenesis, translation, and ATP generation being 346 upregulated on the scalp (Figure 3C, Supplementary Data 2). These observations 347 highlight how skin commensals may employ different gene expression strategies to 348 thrive across body sites, including using specialized nutrient sources (M. restricta) and 349 responding to energetic considerations (C. acnes). 350 The transcriptomes of various skin microbes such as S. epidermidis have been 351 well studied in vitro41,42 but not in vivo, raising important questions about how closely 352 these culture models reflect their actual behaviour on human skin . Skin 353 metatranscriptomics enabled the c omparison of S. epidermidis gene expression in 354 vivo versus three in vitro conditions comprising of cultures at log phase, stationary 355 phase and subjected to osmotic stress , using the same RNA extraction and 356 sequencing protocols (Figure 3D). Gene expression profiles between in vivo and in 357 vitro conditions were clearly separated within the same experimental batch, or across 358 different studies (generated in -house, Avican et al 42 and Wang et al 41) after batch 359 correction (Figure 3E left, Supplementary Figure 11 ). Th ere were many more 360 differentially expressed genes between on -skin versus laboratory -grown conditions 361 (n=1108-1664; Figure 3 E right , Supplementary Data 3, Supplementary Data 6; 362 Methods), relative to the modest changes between moist toe webs and sebaceous 363 cheek/scalp sites in vivo (n=64 differentially expressed genes; Supplementary Data 364 4). Gene set enrichment analysis (GSEA) for in vivo vs in vitro comparisons revealed 365 over 300 gene sets which were significantly enriched in either condition, with most 366 differences observed between sebaceous sites and various in vitro conditions 367 (Supplementary Figure 12, Supplementary Data 5; Methods). While S. epidermidis 368 cultures in vitro upregulated genes in pathways primarily related to carbohydrate 369 metabolism, reflecting the abundance of sugars in rich media , enriched pathways in 370 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint sebaceous sites skewed towards sulphur metabolism, peptide and vitamin 371 biosynthesis, reflective of the more complex nutritional landscape in vivo 372 (Supplementary Figure 1 2). In contrast, there were far fewer significantly enriched 373 gene sets/pathways in toe web sites relative to in vitro conditions despite having 374 similar numbers of differentially expressed genes a s sebaceous versus in vitro 375 comparisons (Figure 3 E, Supplementary Data 3, Supplementary Data 6). 376 Nonetheless, sets of genes involved in metal ion homeostasis (GO:0030001) and 377 copper ion binding (GO:0005507) were consistently upregulated in S. epidermidis 378 colonizing toe webs relative to cultures at log phase or exposed to osmotic stress 379 (Supplementary Figure 13A, Supplementary Data 7). These upregulated genes in 380 the toe webs typically include pumps (e.g. efflux systems and P -type ATPases 43) 381 which couple ATP hydrolysis with the influx and efflux of different substrates and metal 382 ions like copper and zinc ( Supplementary Figure 13B). The upregulation of these 383 pumps could be crucial for survival in toe webs, a “closed" environment rich in sweat 384 containing trace metals excreted by the host44,45, where maintaining a fine balance of 385 intracellular metal levels is essential for microbial cells46. 386 To contextualize overall differences in metabolic activity between in vivo and in 387 vitro conditions via integration through metabolic networks, we performed flux balance 388 analysis (FBA) based on genome-scale metabolic models constrained with our 389 metatranscriptomic data47 (Methods). S. epidermidis showed distinct metabolic fluxes 390 under all the conditions tested (83 metabolic reactions differing across all conditions, 391 PERMANOVA Adonis R 2=0.54, p-value<0.001), with a clear distinction observed for 392 the in vitro conditions (log phase, osmotic stress and stationary phase) and even 393 between in vivo conditions (sebaceous and toe web; Supplementary Figure 14A and 394 Supplementary Data 8). For example, fluxes associated with the production and 395 export of propionate, a key short chain fatty acid impacting skin barrier function and 396 immunity48, were higher in vivo for S. epidermidis compared to in vitro conditions 397 (Figure 4A). Analysis of fluxes through glycolysis and pyruvate generation revealed 398 that unlike under in vitro conditions, S. epidermidis predominantly generates the 399 glycolytic intermediate glyceraldehyde-3-phosphate (G3P) via the pentose phosphate 400 pathway (TALA, PFK_3) and the uptake of L -lactate (L_LACD2) to produce pyruvate 401 for energy generation (Figure 4B). Pyruvate production in vivo is further supported by 402 an active cataplerosis reaction (PPCKr), which generates intermediates contributing 403 to pyruvate synthesis (Figure 4B). Under in vivo conditions, the NADH5 reaction which 404 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint regenerates NAD, showed no flux ( Supplementary Figure 14B). Instead, NAD 405 regeneration might be driven by other compensatory cyclic reactions such as those 406 involving 4-aminobutyrate consumption (ABUTR) and produc tion (ABUTD) 407 (Supplementary Figure 14C). While significant differences were identified in 58 408 reactions between the two in vivo conditions, overall C. acnes exhibited lesser 409 variation in metabolic flux between the scalp and cheek (PERMANOVA Adonis 410 R2=0.053, p-value<0.001; Supplementary Figure 15 and Supplementary Data 9). 411 On both scalp and cheeks, C. acnes displayed an active Wood-Werkman cycle which 412 facilitates propionic acid production49, with higher propionic acid production and export 413 observed on the cheek compared to the scalp (GLM, adjusted p-value<0.001; Figure 414 4C). In vivo site-specific adaptations of C. acnes were also observed in the metabolism 415 of specific amino acids . Flux analysis revealed that glutamate, a critical amino acid 416 integral to linking key carbon and nitrogen pathways 50, was produced via distinct 417 metabolic routes in the scalp and cheek. On the scalp, glutamate production primarily 418 occurred through histidine metabolism (Figure 4 D), whereas on the cheek, it was 419 synthesized predominantly through proline metabolism ( Figure 4 D), emphasizing 420 niche-specific adaptation of metabolic pathways based on resource availability. 421 Overall, these results highlight that species-specific analysis of metatranscriptomic 422 data can be invaluable for identifying metabolic requirements for organismal survival 423 on skin that would not otherwise be reflected in typical in vitro culture conditions. 424 Gene level analysis identifie s key antimicrobial functions and interactions in 425 vivo 426 As antimicrobial peptides and proteins can further shape niche adaptation and inter -427 species interactions, we searched for actively expressed genes that play a role in 428 microbial competition in vivo . Antimicrobial peptides (AMPs) and proteins were 429 sensitively detected by scoring “hit” database genes with profile hidden Markov Models 430 (pHMMs; Methods). Skin microbes showed in vivo expression of a diversity of 431 bacteriocins, phenol soluble modulins, enzymes that generate free radicals, and auto-432 inducing peptides (AIPs) , representing classes of antimicrobial products that have 433 been validated using in vitro or ex vivo models7,8 (Figure 5A). We also determined 434 which species in our dataset expressed various classes of antimicrobial products 435 (Figure 5B), motivated by the importance of identifying candidate species or bio -436 actives for antimicrobial therapies and further experimental follow -up. AIPs were 437 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint expressed by different Staphylococcus species, consistent with the role s of these 438 peptides in Staphylococcal quorum sensing and inter-staphylococcal competition8. Of 439 particular interest are a diverse class of peptides called bacteriocins, which are 440 secreted by bacteria or archaea to inhibit the growth or activities of other microbes in 441 the community 51. For example, some individuals were colonized by strains of S. 442 hominis and S. epidermidis which expressed two peptides of the lacticin 481 family 443 that were recently characterized to have broad antimicrobial activity against a range 444 of gram -positive bacteria 52 (Figure 5B). In contrast, relatively few individuals 445 harboured Staphylococcal species with measurable expression of the well -446 characterized gallidermin/nisin family of bacteriocins which c ould be due to strain 447 differences or its documented autotoxicity53. 448 There were notable examples of site -specific distribution of antimicrobial 449 products. For example, transcripts belonging to cyclic lactone AIPs, thiazolylpeptides 450 and their precursors were more frequently detected on sebaceous sites like the scalp 451 and cheeks compared to other sites (Fisher’s exact test, adjusted p -value<0.001; 452 Figure 5A ). In some cases, site -specific expression was detected despite similar 453 abundances of their host genomes in the community. For example, expression of the 454 lacticin 481 family of lantibiotics by S. hominis and S. epidermidis was not evenly 455 distributed across skin sites, with higher frequencies at the volar forearms and toe 456 webs (Fisher’s exact test, adjusted p-value<0.05; Figure 5A). However, there were no 457 statistically significant differences in metagenomic or total RNA abundances of these 458 two species between lacticin 481 expressors versus non-expressors (Wilcoxon rank 459 sum test, p-value>0.05; Supplementary Figure 16A-B). Such expression variability 460 could be due to strain or environmentally driven differences. In contrast, thiopeptide 461 expression on skin was associated with increased metagenomic and total RNA 462 abundances of C. acnes, pointing to this species as an important source of thiopeptide 463 AMPs (Wilcoxon rank sum test, p-value<0.05; Supplementary Figure 16C-D). 464 Several of the skin microbial bacteriocins detected in this study remain 465 uncharacterized, representing an untapped source for development of new 466 antimicrobials. Thiopeptides are important bio -actives that can shape skin microbial 467 communities, e.g. cutimycin is secreted by C. acnes colonizing the hair follicles and 468 can inhibit the growth of Staphylococcus species2. Two different thiopeptides 469 (MET_03151623 and MET_02967399) were expressed by C. acnes in vivo 470 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint (Supplementary Figure 17). While the thiopeptide MET_03151623 was identical to 471 cutimycin, the other thiopeptide (MET_02967399) represents a putative bacteriocin 472 with uncharacterized properties because it shared <50% amino acid identity with 473 cutimycin and other members of the thiopeptide family (Supplementary Figure 17). 474 Some Cutibacterium and Corynebacterium strains also expressed peptides of the 475 lactococcin 972 family, but most of them had primary sequences which differed from 476 those characterized in bacilli such as Lactococcus (<40% identity ; Figure 5B, 477 Supplementary Figure 18). Unexpectedly, several individuals had strains of S. 478 epidermidis and Staphylococcus pettenkoferi colonizing the toe webs which 479 expressed putative bacteriocins homologous to those of the halocin family (Figure 5B, 480 Supplementary Figure 19). Previously, members of this family were only described 481 in halophilic archaea and bacteria 54,55. These examples highlight the unexplored 482 landscape of antimicrobial mechanisms deployed by various skin microbial species 483 and strains to thrive in their in vivo niches. 484 Finally, we leveraged paired metatranscriptomics and metagenomics to identify 485 putative host -microbe and microbe-microbe interactions . Despite the relatively low 486 fraction of human reads in our metatranscriptomic libraries ( median <3%), a median 487 of 0.8 million human reads per library was obtained, similar to the depth of sequencing 488 commonly seen in single -cell transcriptomic s studies 56,57 (Supplementary Figure 489 20A). Assessing our swab-based metatranscriptomics libraries relative to previous 490 biopsy-based studies 58,59 highlighted that a substantial proportion of human reads 491 were exonic (median 31% versus <15% in a previous biopsy-based study59) and the 492 proportion of intergenic reads was low (median 18%) relative to transcriptomic data 493 generated in-house from skin biopsies (median 32%; Supplementary Figure 20B). 494 Gene Set Variation Analysis60 (GSVA) was used to estimate sample-specific immune 495 pathway activities to find associations between host immunity markers and 496 metagenomic abundances of different microbes. This analysis revealed three 497 statistically significant relationships, all of which involved Staphylococcus capitis 498 (FDR<0.25; Supplementary Data 10; Methods). The IL -6/JAK/STAT3 and Toll -like 499 receptor signalling pathways which are involved in Th17 function61 and microbial 500 detection by the immune system respectively, were positively correlated with 501 Staphylococcus capitis abundances on cheeks ( Supplementary Figure 21). This is 502 consistent with observations showing that high Staphylococcus capitis abundances on 503 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint skin are associated with an IL-17 dominated immune profile in skin disease62. Testing 504 this hypothesis further, we noted that exposure to the supernatants of different strains 505 of S. capitis led to consistently higher levels of both pro-IL-1B and cleaved IL -1B in 506 human keratinocytes relative to th e induction by other commensal Staphylococcal 507 species, in line with our findings that S. capitis is associated with induction of specific 508 immune responses on skin (post-hoc Dunn’s test, adjusted p-value<0.05; Figure 5D). 509 To identify putative m icrobe-microbe interactions , the abundances of 510 transcripts encoding proteins entering the secretory pathway in one species were 511 correlated with organismal abundances of another species within the same skin site 512 (Methods). A total of 13 different combinations of common skin commensals were 513 tested transcriptome-wide and >30 significant associations were identified (FDR<0.1; 514 Supplementary Data 10). Notably, transcript abundances of a protein (DNF11_2196) 515 predicted to enter the secretory pathway in Malassezia restricta were strongly 516 negatively correlated with normalized C. acnes abundances on scalp (Spearman’s ρ<-517 0.7, adjusted p -value<0.05; Supplementary Figure 22A). DNF11_2196 is a poorly 518 characterized gene (https://alphafold.ebi.ac.uk/entry/A0A3G2S5R5) with low primary 519 sequence identity to the top hit in the PDB100 structure database (28.9% identity, 81.2% 520 query sequence coverage ; Foldseek web server, 3Di/AA mode). This makes it 521 challenging to infer function based on sequence-based homology alone. However, the 522 availability of accurate protein structure folding and searching algorithms enables 523 functional inferences based on similarities in 3D structure 63,64. There was greater 524 similarity at the protein structural level indicated by both Foldseek (TM-Score=0.77894, 525 RMSD=2.36) and Dali (Z -score=14.5) between DNF11_2196 and the structure of a 526 Streptomyces papain inhibitor (5ntb-B), which has antimicrobial properties due to its 527 inhibition of bacterial cysteine proteases 65 (Figure 5 C). Importantly, the negative 528 correlation between DNF11_2196 and C. acnes abundances could not be explained 529 by an inverse relationship between the organismal abundances of M. restricta and C. 530 acnes (Supplementary Figure 22B-D). A similar analysis done in cheek samples 531 showed that Cutibacterium granulosum abundances were positively correlated with 532 expression levels of a C. acnes triacylglycerol lipase, independent of relationships 533 between organismal abundances , indicating a potential symbiotic interaction and 534 mechanism (Spearman’s ρ>0.7, adjusted p-value<0.05; Supplementary Figure 22E-535 G). Overall, these results showcase that skin metatranscriptomic and metagenomic 536 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint data can be integrated to identify sets of expressed genes mediating host-microbe and 537 microbe-microbe interactions in vivo, and thus help to prioritize candidate bioactive 538 molecules or signalling pathways for further characterization and validation in model 539 systems and clinical trials. 540

541 We report here the development of a robust experimental and computational workflow 542 tailored for metatranscriptomics of diverse skin sites, circumventing the challenges of 543 relatively low biomass by optimising the sampling, extraction, rRNA depletion, 544 contamination removal , and functional classification steps (Figure 1A -C; 545 Supplementary Figure 1-2). By applying this to a cohort of healthy individuals across 546 multiple skin sites, we reveal the landscape of active species, microbial signatures of 547 adaptation to their niches, as well as several key microbe-microbe and host-microbe 548 interactions in vivo. We note that unlike human stool metatranscriptomes 13 and skin 549 metagenomes, skin metatranscriptomes of the stratum corneum are predominantly 550 composed of microbial reads (>97% of all reads before de-duplication). Consequently, 551 skin metatranscriptomes sequenced to relatively modest sequencing depths (5-10 552 million paired-end reads) can capture most active microbial functions at the community 553 level based on rarefaction analysis (Supplementary Figure 3 ). The abundance of 554 microbial reads means that skin metatranscriptomics can be feasibly deployed for 555 population scale or longitudinal studies, especially as sequencing costs continue to 556 decrease with newer platforms. Our protocol is compatible with a diversity of skin sites 557 across individuals and non-invasive sampling using commercially available swabs 558 makes deployment suitable for clinics to investigate disease cohorts. The moderate-559 to-high success rate of library construction from limited input, along with strong 560 technical reproducibility and even coverage of mapped metatranscriptomic reads 561 across microbial gene bodies, enables robust in vivo gene expression measurements 562 and the assessment of changes in microbial activity in highly heterogeneous cohorts 563 or over time. Given the very low biomass in the antecubital fossa and volar forearm 564 relative to other skin sites, further experimental improvements could be explored to 565 decrease the loss of nucleic acids during the extraction process , such as the use of 566 depletable carrier RNAs and single -cell RNA-seq library kits that can tolerate lower 567 inputs and hence increase robustness. Collapsing of read pairs with identical 5’ and 3’ 568 ends is a conservative approach to de -duplication which removes a substantial 569 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint proportion of microbial reads. More precise de -duplication methods based on unique 570 molecular identifiers (UMIs) can be incorporated into this workflow to better distinguish 571 PCR and biological duplicates66. While we optimized this workflow using swabs which 572 sample the stratum corneum for their ease of use, o ur extraction and computational 573 workflows are also suitable for profiling metatranscriptomes from deeper skin layers 574 obtained from other sampling modalities such as biopsies, pore strips 16 or follicular 575 extracts67. Such studies will be necessary to enrich for signatures of immune cells 576 which typically reside beneath the stratum corneum 68, or for functional analysis of 577 different Cutibacterium acnes strains which have been shown to colonize individual 578 pores in a clonal fashion despite their co -existence across the skin surface 69. 579 Altogether, our workflow is a widely-applicable and clinically tractable approach to 580 survey microbial activities on skin surfaces in vivo. 581 A key observation in our work is that skin metatranscriptomes yield distinct 582 information about microbial activities compared to their matched metagenomes , 583 similar to previously published studies in other contexts such as from human stool13 584 or ocean water 14. Surprisingly, core microbial pathways on skin were expressed by 585 relatively few species (low metatranscriptomic alpha diversity) in the community 586 (Figure 2C ), in contrast to stool metatranscriptomes which harboured numerous 587 “housekeeping” pathways which were expressed by most species in the community13. 588 This implies that unlike stool communities, a comparatively smaller fraction of 589 microbes on skin surfaces are active, with the remainder being either quiescent or 590 dead. Akin to our observations on skin, l ow metatranscriptomic alpha diversity of 591 microbial pathways has also been previously described in nasal cavities and vaginal 592 surfaces9, possibly reflecting nutrient scarcity in these environments which may 593 support the robust growth of only a limited number of species. Of note, Malassezia 594 and Staphyloccocus species had outsized contributions to skin metatranscriptomes, 595 with disproportionally higher RNA abundances relative to their DNA. This can be 596 attributed to a combination of factors such as cell size and bioactivity. As eukaryotic 597 fungi, Malassezia cells have a volume and biomass that are at least two orders of 598 magnitude greater than the average bacterial cell 28. This means that metagenomes, 599 which measure genome copy numbers, likely underestimate the contribution of 600 Malassezia spp. to the functional potential of the skin microbiome. Malassezia species 601 are important sources of secreted lipases, proteases and metabolites that can shape 602 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint host and microbial activities on skin7,70. For example, Malassezia phospholipases can 603 breakdown host lipids to generate polyunsaturated fatty acids (PUFAs) such as 604 arachidonic acid which are potent inflammatory mediators71. The relatively high activity 605 of various Staphylococcus species was unexpected, as they usually make up only a 606 small proportion of bacterial communities on skin, except in sites such as toe webs. 607 This could be reflective of the metabolic versatility of Staphylococcus spp. for survival 608 in diverse skin microenvironments due to their ability to utilize a range of carbon 609 sources, amino acids and lipids for their energetic and cellular demands 72–74. Certain 610 staphylococci such as S. epidermidis have also evolved a range of strategies to 611 escape suppression by the host immune system such as intracellular localization75 or 612 by expressing antigens that lead to commensal -specific T cell responses 76. Our 613 workflow will thus be a useful tool to study how Staphylococcus species contribute to 614 skin health, especially given reported associations of these microbes with skin 615 phenotypes such as malodor77, itch6 and eczema10. 616 Our work represents the first attempt to leverage metatranscriptomics to 617 characterize the metabolic pathways required by different skin microbes to thrive in 618 their in vivo niches. Besides identifying significantly enriched metabolic pathways 619 computed from differential expression analysis, we further extended this approach by 620 integrating transcriptomic data with flux balance analysis (FBA)78. The latter approach 621 models organismal activity as a system of all predicted metabolic reactions instead of 622 examining individual pathways in isolation , allowing us to identify metabolic 623 requirements or dependencies necessary for maximizing biomass in vivo which cannot 624 be readily determined from in vitro co-culture models 79. Interestingly, while the 625 lipophilic fungi Malassezia restricta was present in the metagenomes of both facial 626 skin (cheeks) and scalp, metatranscriptomics revealed differential enrichment of 627 fungal pathways that metabolize distinct classes of host lipids between the two sites . 628 This observation highlights the importance of using metatranscriptomics to study 629 organismal activities or phenotypes, even when there are no significant differences in 630 organismal abundances. When applied to microbes growing under different in vivo 631 versus in vitro conditions, transcriptome aware FBA revealed that the pentose 632 phosphate pathway (PPP) and lactate metabolism were preferentially utilized by S. 633 epidermidis to generate pyruvate as a respiratory substrate in vivo in contrast to in 634 vitro conditions. This observation is consistent with previous reports showing that the 635 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint PPP is crucial for energy production, biofilm formation and virulence in staphylococci80. 636 Interestingly, our results also suggest that besides Cutibacterium acnes, commensal 637 S. epidermidis may be an other in vivo exporter of the short -chain fatty acid (SCFA) 638 propionate, which has been linked to immunomodulation in keratinocyte and sebocyte 639 cell lines81. Further studies are needed to determine if propionate mediates crosstalk 640 between skin commensals and host immune cells in vivo and whether enhanced 641 propionate export by S. epidermidis is a consequence of glucose limitation and/or 642 cysteine/methionine bioavailability ( Figure 4A ). Differing fluxes through metabolic 643 pathways for energy generation and metabolite export in the same organism under in 644 vivo versus in vitro conditions can aid in designing strategies to coax microbes to 645 export desirable metabolites such as SCFAs. This also highlights the need for model 646 systems that better recapitulate the cutaneous microenvironment30 when investigating 647 mechanisms by which commensals such as S. epidermidis contribute to skin health 648 and disease. 649 We are the first study to generate in vivo gene expression data across a range 650 of microbial species, individuals and skin sites to identify key antimicrobial functions 651 and interactions . Although there are reports of skin commensals secreting 652 antimicrobial compounds such as peptides or enzymes which can inhibit pathogens 653 like Staphylococcus aureus, these studies were mostly conducted with in vitro or ex 654 vivo models7,8. We confirmed that many of these products , such as the thiopeptide 655 cutimycin2 and peptides of the lacticin 481 family82 which were described in previous 656 studies, were also expressed in vivo on some individuals by Cutibacterium acnes and 657 commensal Staphylococcus species respectively. Importantly, our dataset shows that 658 the skin microbiome is a rich source of additional antimicrobial proteins/peptides 659 whose activities and specificities remain uncharacterized. One intriguing example is a 660 class of proteins harbouring a C terminal cysteine-rich region that is homologous to a 661 class of bacteriocins known as the halocins. This class of bacteriocins have largely 662 been characterized only in archaea but hypothetical proteins harbouring homologous 663 regions have been found in the genomes of many Staphylococcus strains 664 (https://www.ncbi.nlm.nih.gov/Structure/cdd/TIGR04449). Our data indicates that 665 genes encoding these proteins with halocin C8 -like domains are expressed by 666 commensal Staphylococcus species in vivo (toe webs), providing a basis for further 667 experimental validation of their antimicrobial properties. Since our workflow captures 668 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint the expression of both host and microbial genes/pathways, correlating those with 669 microbial DNA abundances and integrating these results with other functional 670 annotations is a promising approach to identifying novel candidates mediating 671 microbe-microbe and host-microbe interactions in vivo. We identified the Malassezia 672 restricta gene DNF11_2196 as one such candidate owing to the strong negative 673 correlation between transcript levels and C. acnes abundance, its predicted entry into 674 the secretory pathway and its structural similarity to known protease inhibitors (Figure 675 5). Further application of this approach to population-scale cohorts could be useful for 676 robust identification of physiologically relevant interactions between skin microbes and 677 with host cells. 678 In conclusion, we have developed a systematic workflow to sample and analyze 679 skin meta transcriptomes and have shown its utili ty for capturing in vivo microbial 680 activities distinct from just DNA abundances, as well as identifying interactions that 681 can shape microbial communities and host responses. Our data serves as a baseline 682 for healthy individuals to enable future comparisons with datasets for disease states 683 and gain deeper insights into microbial as well as host pathways that could be 684 leveraged for diagnostic and therapeutic applications. 685 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint

686 Subject recruitment 687 The study was conducted in the Genome Institute of Singapore (GIS). All associated 688 protocols for this study were approved by the Agency for Science, Technology and 689 Research Institutional Review Board ( A*STAR IRB reference number 2021-094) on 690 September 8th, 2021, and renewed until September 7th, 2025. All subjects recruited 691 in this study were of Singaporean nationality or permanent residents, aged 21-65 and 692 reported no skin disease at the time of sampling . Five subjects were recruited for 693 cross-sectional and three for longitudinal analysis in a pilot cohort, respectively. For 694 longitudinal analysis, subjects were sampled for 3 visits over 3 consecutive days with 695 a minimum interval of 24 hours. For cross -sectional analysis, 27 subjects were 696 recruited for the full cohort. All subjects were required to abstain from showering for at 697 least 12 hours before sampling. Human skin biops ies from A*STAR Skin Research 698 Laboratories (A*SRL) of 6 and 8mm were obtained from healthy skin donors, with 699 approval from the NHG domain specific review board (NHG DSRB 2017/00224, NHG 700 DSRB 2018/00945) at National University Health System (NUHS) and NSC domain 701 specific review board (NSC DSRB 2019/00806) at National Skin Centre. 702 Optimization of RNA Extraction with pilot cohort 703 Different sample collection methods, bead tubes, purification methods and various 704 combinations of each were tested and are detailed in Supplementary Data 1. The 705 different sample collection modes tested were with FLOQSwabs® (Copan Diagnostics; 706 Cat. No. 502CS01) and tape discs (Cuderm; D-Squame D100). For sample collection 707 with tape discs, tape discs were disinfected with 70% Ethanol for at least 1 minute and 708 dried. Tape discs were applied to skin sites, peeled off and reapplied for a total of 50 709 iterations. Tape discs were then inserted into bead tubes, with the sticky surface facing 710 the inside of the tube. For sample collection with swabs, swabs were first wetted with 711 1× phosphate buffer saline (PBS). The moistened swab was rotated and rubbed with 712 constant pressure in a zig-zag pattern. This procedure was also repeated at an angle 713 of 90° to the first rub, for a total of 1 min ute. Swab contents were either dislodged by 714 vigorous stirring in 800µL or 1300µL of DNA/RNA Shield ( Zymo Research; Cat. No. 715 ZYR.R1100) or were used directly for extraction. RNA extraction with RNeasy 716 PowerMicrobiome kit ( Qiagen; Cat. No. 26000) was carried out according to 717 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint manufacturer’s instructions. RNA extraction with the Quick DNA/RNA Microprep Kit 718 (Zymo Research; Cat. No. D7005) or Norgen BioTek RNA/DNA Purification Micro Kit 719 (Norgen BioTek; Cat. No. 48700) w as carried out according to their respective 720 manufacturer’s instructions. Skin swabs were dislodged into DNA/RNA shield and 721 bead-beating was carried out with ZR BashingBead Lysis tubes (0.1mm & 0.5mm), 722 with 3 cycles of 1 min at 6m/s with 5 min interval on ice on the FastPrep -24 (MP 723 Biomedicals) homogenizer. Lysates were used as inputs into these kits. 724 Samples which were extracted with hot sodium dodecyl sulfate (SDS) and hot 725 acid phenol-chloroform were processed as follows. Skin samples were dislodged from 726 swabs into DNA/RNA shield (2× concentrate; Zymo Research; Cat. No. R1200) and 727 suspensions were added to 0.5× volume of SDS lysis solution (2% SDS, 16mM EDTA) 728 pre-heated to 100°C. These mixtures were incubated at 100°C for 5 min with periodic 729 mixing. Lysates were added to 1× volume of Acid-Phenol:Chloroform:Isoamyl alcohol, 730 pH 4.5 (125:24:1) which was pre-heated to 65°C (Invitrogen; Cat. No. AM9720). This 731 was mixed by vortexing and further incubated at 65°C for 10 min on a thermomixer 732 with constant shaking at 1400rpm. Samples were transferred to 2mL screw cap tubes 733 containing 550mg of autoclaved and air -dried 0.5mm zirconia/silica beads ( BioSpec 734 Products Inc.; Cat. No. 11079105Z). Bead-beating was carried out with 3 cycles of 1 735 min at 6m/s with 5 min interval on ice, using the FastPrep -24 homogeni zer (MP 736 Biomedicals). Samples were centrifuged for 5 min at 16,000g at 4°C and the aqueous 737 phase was transferred to a new tube. 0.1× volume of 3M sodium acetate (pH5.5) 738 (Invitrogen; Cat. No. AM9740), 2× volume of 100% ethanol and 1µL of linear 739 acrylamide (5mg/ml; Invitrogen; Cat. No. AM9520) were added to the aqueous phase. 740 RNA was allowed to precipitate overnight at -20°C and then centrifuged at 16,000g at 741 4°C for 5 min to obtain the pellet. RNA pellets were washed with 1mL of 70% ethanol 742 and air dried for 5 min. Purified RNA was resuspended in RNase -free water, treated 743 with an RNase -Free DNase Set ( Qiagen; Cat. No. 79254) and cleaned up using 744 RNeasy Min-Elute Clean up Kit (Qiagen; Cat. No. 74204). 745 Samples which were extracted with hot phenol were processed as follows. Skin 746 samples were dislodged from swabs into 500µL DNA/RNA shield ( Zymo Research; 747 Cat. No. R1100). A volume of 800µL of phenol solution (Sigma Aldritch; Cat. No. P4682) 748 pre-warmed to 65°C was added to the samples and gently mixed by pipetting without 749 creating bubbles. The phenol-sample mixture was transferred to a tube containing 750 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint 800µL of sodium acetate buffer solution (SAB ; 50mM sodium acetate, 10mM EDTA, 751 pH 5.2). The pre-lysate was transferred into a 2mL screw cap tube containing 550mg 752 of 0.5mm zirconia/silica beads ( BioSpec Products, Inc.; Cat. No. 11079105Z). Bead-753 beating on the FastPrep-24 (MP Biomedicals) homogenizer was carried out at 3 cycles 754 of 1 min at 6.0m/s with 5 min intervals on ice. Lysates were incubated at 65°C for 30 755 min on a thermomixer and mixed at 1400rpm for 10 s every 1 min. Samples were 756 centrifuged at max imum speed for 10 min at 4°C. The aqueous layer was added to 757 750µL of pre-warmed (65°C) phenol and incubated at 65°C for 5 min on a thermomixer, 758 with mixing at 1400rpm for 10 s every 1 min. Samples were centrifuged at 16,000g for 759 10 min at 4°C and the aqueous layer was transferred to a tube containing 750µL of 760 Acid-Phenol:Chloroform:Isoamyl alcohol , pH 4.5 (with IAA, 125:24:1) which was 761 prewarmed to 65°C (Invitrogen; Cat. No. AM9720). The mixture was vortexed for 30 s 762 and left at room temperature for 1 min. Samples were centrifuged at maximum speed 763 for 10 min at 4°C. The aqueous layer was transferred to a new tube and 70µL of 3M 764 sodium acetate and 700µL of isopropanol was added. The mixture was incubated at 765 room temperature for 30 min and centrifuged at 16,000g for 30 min at 4°C to pellet the 766 RNA. The supernatant was removed, and the RNA pellet was washed with 70% cold 767 ethanol. Samples were centrifuged at 16,000g for 5 min, the ethanol was aspirated off 768 and the pellet was air dried. The RNA pellet was then dissolved in 20µL of nuclease 769 free water, treated with RNase-Free DNase Set (Qiagen; Cat. No. 79254) and cleaned 770 up using RNeasy Min-Elute Clean up Kit (Qiagen; Cat. No. 74204). 771 Samples were extracted with the Direct-zol™ RNA MicroPrep kit (Zymo 772 Research, Cat. No. ZYR.R2063) with Lysing matrix E ( MP Biomedicals ; Cat. No. 773 116914500) or ZR BashingBead Lysis tubes (0.1mm & 0.5mm; Zymo Research; Cat. 774 No. S6012) or zirconia/silica beads ( BioSpec Products; Cat. No. 11079105Z). Tape 775 discs were directly added into bead tubes, while skin swabs or bacteria l cells were 776 dislodged into DNA/RNA shield. RNA purification was carried out according to the 777 manufacturer’s protocol, except for some modifications such as the addition of 0.05% 778 Tween 20 and using either 3 or 4 cycles of bead -beating, as detailed in 779 Supplementary Data 1. Purified RNA was treated with RNase -Free DNase Set 780 (Qiagen; Cat. No. 79254) and cleaned up using the RNeasy Min -Elute Clean up Kit 781 (Qiagen; Cat. No. 74204). 782 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint Sample collection for full cohort 783 Skin s amples were collected using FLOQSwabs® (Copan Diagnostics ; Cat. No. 784 502CS01) from five different skin sites (scalp, cheek, antecubital fossa, volar forearm, 785 and toe web) from each subject (n=135 metagenomes and n=135 786 metatranscriptomes). For each skin site, three swabs were collected from the left and 787 right side of the body and combined in a tube, except for the scalp, where only three 788 swabs were used in total. Each swab was submerged in 1× phosphate buffer saline 789 (PBS), and the excess solution was removed by pressing the swab against the tube 790 wall. For each skin site, the moistened swab was rotated and rubbed, with constant 791 pressure applied, in a zig-zag pattern and was repeated at an angle of 90° to the first 792 rub, for a total of 1 min. The contents of the swab were dislodged by stirring vigorously 793 in either 800µL or 1300µL of DNA/RNA Shield (Zymo Research; Cat. No. ZYR.R1100) 794 for scalp and other skin sites respectively. Swabs were submerged in DNA/RNA shield 795 for 5 min at room temperature and stirred vigorously again. Excess solution on each 796 swab was collected by pressing the swab against the wall of the tube. This was 797 repeated for the remaining swabs for each site. Negative controls (n=7) for each batch 798 were collected by dislodging three swabs in DNA/RNA Shield (Zymo Research ; Cat. 799 No. ZYR.R1100) without sampling the skin. Each sample or control was split into 2 800 portions for RNA (approximately 600µL) and DNA (approximately 200µL) extraction. 801 All samples were stored at -80°C prior to nucleic acid extraction. 802 RNA extraction for metatranscriptomics (Direct-zol method) 803 RNA was extracted using the Direct -zol™ RNA MicroPrep kit (Zymo Research ; Cat. 804 No. ZYR.R2063). Compared to other extraction methods and kits, this approach was 805 found to be the best performing in terms of RNA yield and RNA integrity 806 (Supplementary Data 1). Bead-beating of samples in TRI Reagent® (Zymo Research; 807 Cat. No. R2050-1-200) was done using ZR BashingBead Lysis Tubes (0.5 & 0.1mm ; 808 Zymo Research; Cat. No. S6012) and Fastprep -24 Instrument (MP Biomedicals) at 809 6.0m/s for a total of 3 min, in 1 min intervals, with 5 min incubation on ice between 810 each interval. Samples were DNase treated and purified according to kit instructions. 811 An additional DNase treatment was carried out by adding 2.5µL of DNase I and 10µL 812 of RDD buffer from RNase -Free DNase Set (Qiagen ; Cat. No. 79254), and 1µL of 813 Recombinant RNasin RNase Inhibitor (10,000U; Promega; Cat. No. N2515) in a total 814 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint volume of 100µL. The mixture was incubated at 37°C for 30 min and purified with 815 RNeasy MinElute Cleanup Kit (Qiagen; Cat. No. 74204) in an elution volume of 14µL 816 of RNase-free water. High Sensitivity RNA ScreenTape analysis (Agilent Technologies; 817 Cat. No. 5067-5579, 5067-5580) was used to assess the quality of RNA and extracted 818 RNAs were stored at -80°C. 819 Mock communities 820 For metagenomic spike-ins, a mock community of 3 different bacteria (Vibrio vulnificus 821 ATCC® 29307 ™, Plesiomonas shigelloides ATCC® 51903 ™, and Listeria 822 monocytogenes ATCC® 35152™) was created by culturing each bacterial species in 823 their respective growth media as recommended by ATCC. A growth curve was 824 established for each bacterial species and the optical density (OD) of each culture was 825 used to estimate the number of colony forming units (CFU s). Equal number of CFUs 826 for each species were resuspended in DNA/RNA shield (Zymo Research ; Cat. No. 827 ZYR.R1100) and pooled together to create a mock community stock of concentration 828 1×105 CFU/µL. Mock community stocks were aliquoted and stored at -80°C. Other 829 mock communities for quality control testing were prepared to consist of 10 different 830 bacterial species commonly part of human skin and gut microbiomes (Supplementary 831 Table 2). Each species was cultured according to the recommended growth conditions 832 and was counted using a hemocytometer. Cultures were centrifuged and each cell 833 pellet was re-suspended in 100-500µL of DNA/RNA Shield (Zymo Research; R1100). 834 Samples were stored at -80°C prior to combining the cells according to the 835 abundances summarized in Supplementary Table 2. 836 DNA extraction for metagenomics (EZ1 method) 837 This approach (EZ1 method) was used for DNA extraction due to relatively poor DNA 838 yields from the Direct-zol method, which is optimized for RNA extraction. Spike-ins of 839 the metagenomic mock community (1.5×104 CFUs) were introduced to each sample 840 prior to DNA extraction. Lysis of samples was carried out by adding 500µL of ATL 841 Buffer (Qiagen; Cat. No. 19076) to the sample and homogenisation in Lysing Matrix E 842 tubes (MP Biomedicals ; Cat. No. 116914500) with a FastPrep-24 Instrument at a 843 speed of 6.0m/s for 40 s, done twice in total. Cell debris was pelleted at 16,000g for 5 844 min and the supernatant was treated with 12µL of Proteinase K at 56°C for 15 min 845 prior to purification with EZ1 DNA Tissue Kit (Qiagen; Cat. No. 953034) using the EZ1 846 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint Advanced XL machine (Qiagen). A Qubit fluorimeter was used to quantify the amount 847 of DNA. 848 Extraction quality control 849 To assess the lysis efficiencies of different protocols, each mock community 850 (Supplementary Table 2) was extracted with EZ1 DNA Tissue Kit (Qiagen ; 953034; 851 EZ1 method) and a separate aliquot with TRI Reagent® (Zymo Research; Cat. No. 852 R2050-1-200; Direct-zol method). DNA was obtained from the Direct-zol method using 853 additional steps after lysis according to the TRIzol™ Reagent (DNA isolation) User 854 Guide. DNA was cleaned up and concentrated with 2× AMPure XP beads (Beckman 855 Coulter; A63882). DNA was eluted in 30µL of Buffer EB (Qiagen; 19086). Qubit RNA 856 High Sensitivity assay kit (Thermo Fisher Scientific; Q32852) was used to check for 857 RNA contamination. Samples extracted by the Direct-Zol method that had detectable 858 levels of RNA were treated with RNase A (100mg/mL ; Qiagen; 19101), cleaned up 859 with 2× AMPure XP beads, and eluted in 30µL of Buffer EB. DNA was used for library 860 preparation and sequenced on a Illumina Novaseq X (2×150 bp reads). This data was 861 used to confirm that metagenomic relative abundances obtained based on DNA 862 extracted from both methods were highly co nsistent across the different mock 863 communities that were tested (Pearson’s R=0.95; Supplementary Figure 23). 864 Preparation of ribosomal RNA depletion mix 865 Five different oligo probe pools (desalted) were ordered from Integrated DNA 866 Technologies (IDT) targeting either Malassezia rRNA or other fungal rRNAs 867 (Supplementary Table 1). The probe pools for Malassezia rRNAs were mixed to form 868 a pool of 868 probes, while probe pools for other fung al rRNAs formed a pool of 630 869 probes. Both oligo pools were resuspended to a concentration of 2 µM/probe. The two 870 oligo probe pools were then mixed together (i.e., Malassezia:other fungi – 4:1) to form 871 a pan -fungal rRNA depletion probe pool. The final rRNA depletion probe pool was 872 made by combining the NEBNext rRNA Depletion Solution from NEBNext rRNA 873 Depletion Kit V2 (Human/Mouse/Rat ; New England Biolabs ; Cat. No. E7405) and 874 NEBNext rRNA Depletion Kit V2 (Bacteria ; New England Biolabs ; Cat. No. E7850) 875 with the custom pan-fungal rRNA depletion probe pool at a volume ratio of 40:9:1, with 876 2 µL of this final probe mix used per sample for rRNA depletion. 877 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint RNA library preparation 878 Human and microbial ribosomal RNAs were depleted from 5–10 ng of total RNA or the 879 entire volume of eluted total RNA . Libraries were prepared according to 880 manufacturer’s instructions using NEBNext Ultra II Directional RNA Library Prep Kit 881 for Illumina (New England Biolabs ; Cat. No. E7760). Depending on the RIN value of 882 RNA, either Section 2 or 3 of the protocol was used. Library enrichment was carried 883 out using NEBNext Multiplex Oligos for Illumina (96 Unique Dual Index Primer Pairs ; 884 New England Biolabs ; Cat. No. E6440) or NEBNext Multiplex Oligos for Illumina 885 (Unique Dual Index UMI Adaptors RNA Set 1; New England Biolabs; Cat. No. E7416) 886 with 14 cycles of enrichment PCR. The q uality of libraries w as assessed by High 887 Sensitivity D1000 ScreenTape Assay (Agilent Technologies ; Cat. No. 5067 -558). 888 Libraries were pooled in equimolar proportions and sequenced on a Illumina HiSeq X 889 Ten system (~35 million 2×150bp read pairs per library). 890 DNA library preparation 891 NEBNext Ultra II FS DNA Library Prep Kit for Illumina (New England Biolabs; Cat. No. 892 E7805) was used according to manufacturer’s instructions with some modifications. A 893 volume of 26µL of DNA was used as input and subjected to 10 min of fragmentation 894 at 37°C. Fragmented DNA was used for adapter ligation and was cleaned up using 895 0.6× volume of AMPure XP Reagent (Beckman Coulter ; Cat. No. A63882). Adaptor 896 ligated DNA was amplified for 12 cycles using NEBNext Multiplex Oligos for Illumina 897 (96 Unique Dual Index Primer Pairs ; New England Biolabs; Cat. No. E6440) and 898 cleaned up with 0.7× volume of AMPure XP Reagent (Beckman Coulter ; Cat. No. 899 A63882). The final library was eluted in 20µL of EB Buffer (Qiagen ; Cat. No. 19086). 900 Quality of libraries were assessed by the High Sensitivity D1000 ScreenTape Assay 901 (Agilent Technologies ; Cat. No. 5067 -558). Libraries were pooled in equimolar 902 proportions and sequenced on a Illumina HiSeq X Ten system (~25 million 2×150bp 903 read pairs per library). 904 Data pre-processing and quality control 905 Short reads from metagenomic and metatranscriptomic libraries were processed using 906 a Nextflow83 pipeline (https://github.com/Chiamh/meta-omics-nf). Quality control and 907 adapter trimming were done using fastp 84 (v0.22.0) with default settings . 908 Metagenomes were further pre -processed by mapping to the hg38 human reference 909 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint genome using BWA-MEM85 (v0.7.10-r789) and reads that failed to map to hg38 were 910 extracted using samtools 86 (v1.13) with parameters −f12 −F256. Human RNA was 911 removed from metatranscriptomes by mapping to hg38 using STAR87 (2.7.9a). Reads 912 originating from microbial ribosomal RNAs (rRNAs) were computationally removed 913 from metatranscriptomes using bbduk.sh (BBMap v38.93) and a k -mer database for 914 rRNAs. Subsequently, microbial RNA reads were de -duplicated using clumpify.sh 915 (BBMap v38.93) with parameters dedupe=t and optical=f. 916 Taxonomic classification 917 Metagenomic reads were classified using Kraken2 24 (v2.1.2) and Bracken 88 (v2.6.1) 918 with parameters --use-names, --paired and --report-minimizer-data. 919 Metatranscriptomic reads were classified with Kraken2 using the same parameters. A 920 50Gbp Kraken2 database built from Ref Seq bacterial, archaeal, viral, fungal and 921 human (hg38) genomes, as well as plasmid sequences was used. This database also 922 contains additional Malassezia assemblies downloaded from NCBI ( Supplementary 923 Data 11). Samples with at least 10,000 paired reads were retained and reads that 924 were still taxonomically assigned to Homo sapiens were removed. Microbial reads 925 were defined as the sum of reads classified as bacteria (taxid 2), archaea (taxid 2157), 926 virus (taxid 10239), and fungi (taxid 4751). 927 False positive species assignments for metagenomic reads were identified and 928 removed using an approach similar to Breitwieser et al 23. Species assignment of 929 metagenomic reads were considered true positives if there were ≥2000 unique 930 Kraken2 minimizers per 1 million microbial reads or had ≥10 read pairs for the species 931 with ≥10× more unique minimizers than read pairs. False positive species assignments 932 for metatranscriptomic reads were identified and removed using a similar approach, 933 with empirically determined minimizer thresholds. Specifically, RNAs from three 934 biological repeats (MHS445, MHS589 and MHS590), each comprising a mix of non -935 skin bacteria Plesiomonas shigelloides, Vibrio vulnificus and Listeria monocytogenes 936 were extracted, processed, and sequenced in the same way as other 937 metatranscriptomes in this study. Random subsamples of these libraries were added 938 to a skin metatranscriptome (MHS413) such that read pairs from P. shigelloides, V. 939 vulnificus and L. monocytogenes were each present at relative abundances of 940 approximately 0.1%, 1% or 10% of all reads in the metatranscriptome. These 941 metatranscriptome mixtures were subsequently rarefied to either 104 or 106 read pairs 942 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint and reads were classified using Kraken2. For each true positive spike -in species (P. 943 shigelloides, V. vulnificus and L. monocytogenes ), false positive classifications of 944 reads to other species belonging to the genera Listeria, Plesiomonas, Vibrio, 945 Pseudomonas, Salmonella and Klebsiella were identified. The log 10 ratio of unique 946 species-specific minimizers per million microbial reads was plotted against log10 paired 947 counts of each species to determine that 10 4 unique minimizers per million microbial 948 reads could distinguish most true positive species assignments from false positives 949 (Supplementary Figure 2C). Species assignment of metatranscriptomic reads were 950 considered true positives if there were ≥10 4 unique Kraken2 minimizers per 1 million 951 microbial reads or had ≥10 read pairs for the species with ≥10 × more unique 952 minimizers than read pairs. Species contributing <0.1% RNA relative abundance were 953 removed from metatranscriptomes to avoid false positives from incorrect read 954 classification. Species contributing <0.1% DNA relative abundance were also removed 955 from each metagenome unless they were also present in the matched 956 metatranscriptome (Supplementary Figure 2C). 957 Kitome removal 958 Potential reagent and laboratory contamination-associated species (the “kitome”) 959 were identified and removed via a multi -step process ( Supplementary Figure 2A). 960 Swab extraction controls (n=7 negative control s) were sequenced from fresh swabs 961 unexposed to human skin. Species with ≥0.1% relative abundance in any negative 962 control constituted an initial list of 70 candidate genera. Each of these 70 genera were 963 classified as kitome candidates if they were not previously reported on skin or in skin 964 diseases, based on hits to the Disbiome database89, the MicrophenoDB90 or a PubMed 965 literature search (last accessed 24th Nov 2022) using the search terms ([Genus name] 966 AND "Skin" AND "microbiome" AND "human"). Pairwise Pearson correlations between 967 all species in the metagenomes or metatranscriptomes and these kitome candidates 968 were calculated using SparCC91. Any species with a strong positive correlation (r≥0.8) 969 with kitome candidate species were also marked as potential kitome contaminants and 970 removed from downstream analyses. We observed that abundances of previously 971 reported92 environmental and kit contaminants such as species of the genera 972 Achromobacter, Bradyrhizobium, Mycolibacterium, Mycobacterium and 973 Brevundimonas were not strongly positively correlated with skin or oral microbes and 974 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint could hence be easily distinguished from them for removal ( Supplementary Figure 975 2B). 976 Functional classification 977 Metagenomic and metatranscriptomic reads were functionally classified based on a 978 similar strategy to that of HUMAnN3 20, using a custom Nextflow pipeline 979 (https://github.com/Chiamh/meta-omics-nf). Reads were first aligned in single -end 980 mode using Bowtie293 (v2.4.4) in –very-sensitive mode, to the Integrated Human Skin 981 Microbial Gene Catalog 19 (IHSMGC) comprising approximately 10.9 million non -982 redundant genes. A coverage filter of 50% across the length of any given hit pangene 983 sequence was used 20. Reads which failed to align to the IHSMGC were mapped 984 against the Uni Ref9094 database (downloaded 9 th June 2021) using Diamond 95 985 (v2.0.12) with parameters blastx, --id 80 --query-cover 90 and --max-target-seqs 1. For 986 searches against the UniRef90 database, only alignments with ≥80% sequence 987 identity, ≥90% query (read) coverage and ≥50% subject (UniRef90 representative 988 sequence) coverage were considered as hits. These thresholds were the same as the 989 defaults used by HUMAnN3 and were previously shown to increase specificity of 990 alignments to the UniRef90 database with a relatively small reduction in sensitivity20. 991 Pangenes and UniRef90 clusters with valid hits after mapping were annotated 992 and grouped into orthologous groups (OGs) using EggNOG mapper96 (v2.1.6) and the 993 EggNOG 5 .0 database97 with parameters -m diamond and --go_evidence all. This 994 provides annotations for Gene Ontology (GO)98,99 and Kyoto Encyclopaedia of Genes 995 and Genomes100 (KEGG) pathway analysis. Gene-level analysis was done for OGs by 996 summarizing the read counts at the level of bacteria (taxid 2) or fungi (taxid 4751). 997 Rarefaction analysis for bacterial and fungal OGs was conducted in R using the 998 ‘rarecurve’ function from the vegan package (v2.6-6.1). 999 Pathway abundance and contributional diversity analysis 1000 Pathway abundances were computed using HUMAnN3 (v3.8), using a custom 1001 structured KEGG module definition file 1002 (https://github.com/CSB5/skin_metatranscriptome), with each definition retrieved 1003 using the KEGG REST API (e.g. https://rest.kegg.jp/get/M00357). Microbial alpha 1004 (Simpson) and beta (Bray -Curtis dissimilarity) diversity for each KEGG module was 1005 calculated as described in Franzosa et al9. Only modules that were core to a skin site 1006 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint (non-zero counts in >75% of individuals) with largely known microbial provenance (0.5 were considered “high” diversity. 1009 Staphyloccocus epidermidis culture experiments 1010 Overnight cultures of S. epidermidis (ATCC 12228) were diluted into 3 volumetric 1011 flasks of culture at OD 600 0.01 and incubated at 37°C for 6h to log phase (OD 600 0.4-1012 0.5). Each flask was split into 3 tubes (biological triplicates) and subjected to different 1013 stress exposures as described in Avican et al 42 (Supplementary Table 2). A portion 1014 (2×500µL) of each culture was transferred into 2mL tubes after stress exposure and 1015 1.5mL of DNA/RNA Shield (Zymo Research; Cat. No. R1100) was added into each 1016 tube and mixed well. Cell mixtures were then centrifuged at 8000g for 10min at 4°C to 1017 pellet cells. Each cell pellet was re -suspended in 500µL of TRI Reagent (Zymo 1018 Research; Cat. No. R2050-1-200) and each technical duplicate was combined into 1 1019 tube with a final volume of 1mL. Samples were stored at -80°C prior to RNA extraction. 1020 Alignment to species-specific pangenomes 1021 Bacterial species-specific analyses w as done by mapping metagenomic and 1022 metatranscriptomic reads to curated pangenomes101 (https://ngdc.cncb.ac.cn/propan/) 1023 of eight commonly found skin microbes ( Staphylococcus aureus, Staphylococcus 1024 epidermidis, Staphylococcus hominis, Staphyloccocus capitis, Cutibacterium acnes, 1025 Cutibacterium modestum, Corynebacterium tuberculostearicum and Corynebacterium 1026 ureicelerivorans), together with decoy genomes of the non -skin microbes 1027 Achromobacter xylosoxidans, Plesiomonas shigelloides, Vibrio vulnificus and Listeria 1028 monocytogenes. All genes in the eight skin species pangenome s were further 1029 clustered using CD -HIT (v4.8.1) at ≥95% protein sequence identity and mutual 1030 alignment coverage ≥90%, following guidelines by Li et al 19. Reads were aligned to 1031 the de-replicated eight skin species pangenome s in single-end mode using Bowtie2 1032 (v2.4.4) in --very-sensitive mode (similar to HUMAnN3) . A coverage filter of 50% 1033 across the length of any given hit pangene sequence was used. Fungal species -1034 specific analyses were done by pseudo-alignment of metatranscriptomic reads to the 1035

transcriptomes of multiple Malassezia species, with their reference 1036 genomes as decoys to minimize occurrences of non -transcriptomic reads being 1037 erroneously counted due to similarities to the annotated transcriptome, following 1038 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint recommendations for Salmon102 (v1.10.1). Unlike fungal reads, bacterial reads were 1039 mapped in single-end mode to account for the organization of ORFs in polycistronic 1040 mRNAs, similar to the approach adopted by HUMAnN3. Single-ended read coverage 1041 over bacterial coding sequences w as computed using picard (v3.1.1) 1042 CollectRnaSeqMetrics with arguments -STRAND 1043 FIRST_READ_TRANSCRIPTION_STRAND. Paired -end read coverage over fungal 1044 transcripts were computed using picard CollectRnaSeqMetrics with arguments -1045 STRAND SECOND_READ_TRANSCRIPTION_STRAND. 1046 Transcriptional activity analysis 1047 An organism’s transcriptional activity was calculated in a similar way as described in 1048 Abu-Ali et al 13. Metagenomic or metatranscriptomic reads were mapped to the 1049 IHSMGC19 pangene catalogue and the UniRef90 database. Read counts for any given 1050 species (including all sub-species) were then divided by the length of each feature in 1051 kilobases to obtain the number of reads per kilobase (RPK). Reads belonging to a 1052 species that could not be mapped to pangene or UniRef90 features were assumed to 1053 belong to an unknown gene of length 1 kilobase for computing RPK. RPKs for each 1054 metagenomic or metatranscriptomic library were then summed and divided by 10 6 to 1055 obtain a per-sample scaling factor. Metatranscriptomic transcripts per million (TPM) or 1056 metagenomic copies per million (CPM) values were computed by dividing a feature’s 1057 RPK with the per -sample scaling factor. The transcriptional activity of a species was 1058 estimated by summing species-level TPM values (RNA) and dividing them by species-1059 level CPM (DNA) values. 1060 Differential expression analysis 1061 Differential expression analysis was done using DESeq2 103 (v1.36.0). Raw counts 1062 were summarized at the level of bacterial (taxid 2) and fungal (taxid 4751) orthologous 1063 groups for gene -level differential expression analysis. The design formula was ‘~ 1064 subject ID + skin site (Sc, Ch, Ac, Vf or Tw) + assay (RNA or DNA) + skin site:assay’. 1065 This design formula accounts for within-subject dependencies and variations in gene 1066 copy numbers in metagenomes while testing for differences in microbial gene 1067 expression between skin sites. Only features (rows) with median read count ≥10 for 1068 both DNA and RNA were kept. Size factors were estimated separately for the 1069 metagenomic and metatranscriptomic count matrices using the ‘poscounts’ function to 1070 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint account for data sparsity. Differential expression analysis at individual species -level 1071 between two in vivo conditions was similarly done, except that only metatranscriptomic 1072 read counts were used from libraries with ≥200,000 species -specific reads, and the 1073 design formula was ‘~ subject + skin site ’. Differential expression analysis for 1074 Staphylococcus epidermidis comparing in vivo versus in vitro growth conditions was 1075 similarly done, except that the design formula was ‘~ batch + growth condition ’ to 1076 account for experimental batch effects when using in vitro RNA-seq data from different 1077 sources. Batch corrected principal component analysis (PCA) plots were derived from 1078 inputs processed with the removeBatchEffect function from limma104 (v3.60.4). 1079 Gene set enrichment analysis 1080 Gene set enrichment analysis (GSEA) for differentially expressed microbial features 1081 was done using clusterProfiler 105 (v4.4.4) with the arguments : eps=0, 1082 nPermSimple=10000 and seed=TRUE. 1083 Integration of metatranscriptomics data with metabolic models 1084 Genome-scale metabolic models (GSMMs) for Staphylococcus epidermidis ATCC 1085 12228 and Propionibacterium acnes KPA171202 were obtained from the AGORA 1086 database110. These models were constrained using COBRApy 111 based on specific 1087 exchange flux values corresponding to the conditions under which the simulations 1088 were performed. Genes in the metabolic models were mapped based on transcript 1089 levels, calculated as the geometric mean of transcript abundance (TPMs) across 1090 replicates or samples under the same growth condition. The integration of these TPM 1091 values into genome -scale models and the subsequent flux balance analysis (FBA) 1092 was carried out in Python (v3.12) using R IPTiDe112 (v3.4.81). Briefly, R IPTiDe 1093 incorporates gene expression data using reaction parsimony, generating context -1094 specific GSMMs. The context -specific models were simulated using FBA to identify 1095 flux distributions within the organism under each condition separately. Non-metric 1096 multidimensional scaling (NMDS) of the Bray -Curtis distances between flux 1097 distributions was performed using functions from the vegan package (v2.6-4) in R 1098 (v4.3.0) to compare the flux profiles across conditions. Differentially abundant 1099 reactions were identified using a generalized linear model (GLM: reactions ~ group), 1100 with reactions showing adjusted p -values below 0.001 and absolute estimate values 1101 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint greater than 2 being considered significantly different between conditions . These 1102 reactions were plotted using ggplot2 (v3.5.1) in R. 1103 Identification and analysis of antimicrobial genes 1104 Hidden Markov Models (HMMs) of various classes of antimicrobial proteins were 1105 downloaded from NCBI (https://www.ncbi.nlm.nih.gov/protfam; Supplementary Data 1106 12). Microbial pangenes and representative sequences from Uni Ref90 gene clusters 1107 were searched against these HMMs using hmmscan from hmmer 106 (v3.3.2). Hits to 1108 a HMM were only kept if both the “seq” and “best one domain” scores were greater 1109 than or equal to the sequence and domain cutoffs given by NCBI. Antimicrobial genes 1110 were considered “present” in metatranscriptomes if their read counts were ≥5, with 1111 ≥50% coverage over the gene body. Multiple sequence alignment (MSA) of microbial 1112 sequences and representative sequences from HMMs was done using MUSCLE107 1113 (v5.1). Matrices of pairwise percentage identities were computed from the MSAs and 1114 coverage statistics were computed using a custom python script 1115 (https://github.com/CSB5/skin_metatranscriptome). 1116 Identification of microbe-gene correlations 1117 SignalP108 (v6.0) was used to classify microbial proteins (features) that can enter the 1118 secretory pathway. Fungal and bacterial proteins were analysed in “fast” mode with 1119 the options --organism “eukarya” or “other” respectively. Microbial features predicted 1120 to enter the secretory pathway were shortlisted for correlation analysis. For a given 1121 pair of microbes at a skin site, pairwise Spearman correlations were computed 1122 between the variance stabilized counts (RNA) of microbial features computed from 1123 DESeq2 and the centered log-ratio109 (clr) transformed counts (averaged across 1000 1124 Monte Carlo instances) of microbial abundances (DNA) computed from aldex2 1125 (v1.28.1) and the aldex2propr function from propr (v2.1.2). Features were considered 1126 significantly correlated only with Spearman’s ρ≥0.7 and FDR adjusted p-value≤0.05. 1127 Structural similarity searches 1128 Protein Data Bank (PDB) files of selected microbial proteins were downloaded from 1129 the AlphaFold Protein Structure Database ( https://alphafold.ebi.ac.uk/). Structural 1130 similarity searches were done using the Foldseek Search server 1131 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint (https://search.foldseek.com/search) in 3Di/AA mode and the Dali server 1132 (http://ekhidna2.biocenter.helsinki.fi/dali/) in PDB search mode. 1133 Microbial strain isolation 1134 Commensal Staphylococcal strains were isolated from skin swabs for 4 healthy donors 1135 across 4 separate body sites: antecubital fossa, axilla, cheek and scalp. Swabs were 1136 immersed in 2.5mL BHI broth (Oxoid), incubated at 37°C at 210rpm for 2 hours before 1137 being spread onto TSA -SB (Thermo Scientific) and Baird Parker (Oxoid) plates and 1138 incubated at 37°C for 24 hours. Colonies were then picked and grown in 4ml BHI broth 1139 for 16 hours overnight at 230rpm and 37°C. 200μL of the cultures were spun down at 1140 5000rpm for 5 minutes. To extract gDNA the pellet was resuspended in 50 μL 1141 QuickExtract (Lucigen). Bacterial suspensions were then heated to 65°C for 6 minutes, 1142 briefly vortexed, then heated to 98°C for 4 minutes and briefly vortexed again. PCR 1143 reactions were setup with 10μL KAPA SYBR FAST (Merck), 6.4μL nuclease free water 1144 (Promega), 3μL extracted gDNA and 0.6μL primer mix. Primer mixes were specific for 1145 4 specific Staphylococcal species (Supplementary Table 1 ). PCR was run on an 1146 AriaMx Real Time PCR System (Agilent Technologies). PCR products were sent for 1147 Sanger sequencing to confirm product sequence s. The identified Staphylococcal 1148 colonies were inoculated into 4ml BHI broth and grown overnight for 16 hours at 1149 230rpm a nd 37°C. Optical density (OD) was measured with a SpectraMax M5 1150 Microplate Reader (Molecular Devices). Overnight cultures were then grown from OD 1151 0.1 in fresh BHI for 26 hours at 230rpm a nd 37°C. Final OD measurements for all 1152 strains after 26 hours were normalized to OD 2.75. Bacterial cultures were spun down 1153 at 5000rpm for 5 minutes, pellets were discarded and supernatants stored at 4°C. 1154 Human keratinocyte cell culture 1155 A N/TERT keratinocyte cell line was used for ELISA experiments. A genetically 1156 modified N/TERT keratinocyte cell line was utilized for Nano-Glo HiBiT experiments 1157 which contained a 33 nucleotide HiBiT tag directly adjacent to the IL1B start codon. 1158 Keratinocytes were cultured in keratinocyte serum free medium (KSFM; Gibco) 1159 supplemented with bovine pituitary extract (BPE) at a final concentration of 20μg/mL, 1160 epidermal growth factor (EGF) at a final concentration of 0.2ng/mL, calcium chloride 1161 at a final concentration of 300 μM and 1:1000 penicillin -streptomycin (Gibco). 1162 Keratinocytes were seeded onto 96 well cell culture plate (Greiner) at a density of 1163 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint 20,000 cells per well and grown for 24 hours in 100 μL KSFM. After 24 hours KSFM 1164 was removed and the keratinocytes were cultured in 100 μL KSFM with either 1 μM 1165 anisomycin (Merck), 0.2% TritonX (Merck), 5% BHI or 5% Staphylococcal 1166 supernatants. Keratinocytes were cultured overnight for 16 hours in treatment 1167 conditions. 1168 Nano-Glo HiBiT for pro-IL-1B measurements 1169 After 16 hours of treatment of the HiBiT N/TERT keratinocytes, 50μL of the KSFM was 1170 transferred to a white-bottom 96 well plate (Sigma-Aldrich) and mixed with 50μL Nano-1171 Glo HiBiT extracellular reagent. Plates were mixed on an orbital shaker for 1 minute, 1172 incubated for 10 minutes at room temperature before luminescence was measured 1173 with a SpectraMax M5 Microplate Reader. 1174 ELISA for cleaved IL-1B 1175 After 16 hours of treatment o f the N/TERT keratinocytes , 100μL of the KSFM was 1176 transferred to a new 96 well cell culture plate and stored at -80 °C until required. 1177 Human IL -1B ELISA were performed in accordance with the manufacturer's 1178 instructions (FineTest), OD measurements were performed with a Spark Multimode 1179 Microplate Reader (Tecan) at 450nm and corrected against 570nm. 1180 Statistical analysis and visualization 1181 Statistical tests and visualizations were done using R (v4.2.0), ggplot (v3.3.6), 1182 EnhancedVolcano (v1.14.0) and ggpubr (v0.4.0). 1183 Data availability 1184 Shotgun metagenomic and metatranscriptomic sequencing data is available from the 1185 European Nucleotide Archive (ENA – https://www.ebi.ac.uk/ena/browser/home) under 1186 project accession number PRJEB82796. All large datasets, Foldseek webserver 1187 outputs, Dali webserver outputs , the Kraken2 database used for taxonomic 1188 classification, and other gene annotation databases are available on Figshare at 1189 https://figshare.com/projects/Skin_metatranscriptomics_manuscript_2024/202683. 1190 Code availability 1191 Source code for scripts used to analyze the data are available at 1192 https://github.com/CSB5/skin_metatranscriptome. 1193 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint

1194 The authors wish to thank Lim Thiam Chye and Nelson Teo from NUHS, and Steven 1195 Thng and Yew Yik Weng from NSC for collecting patient skin biopsies. This work is 1196 supported by the Asian Skin Microbiome Programme 2.0 (Industry Alignment Fund 1197 Pre-Positioning; H22J1a0040), Agency for Science, Technology and Research BMRC 1198 EDB IAF -PP grants – H17/01/a0/004 and Agency for Science, Technology and 1199 Research BMRC Central Research Funds (ATR) and a National Medical Research 1200 Council (NMRC) Clinician Scientist-Individual Research Grant (CIRG23jul-0018). This 1201 work was also supported by the A*STAR Computational Resource Centre through the 1202 use of its high -performance computing facilities. The N/TERT keratinocyte cell lines 1203 were provided by the Zhong Lab, from LKC-Medicine to A*SRL. 1204 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint Figure Legends 1205 Figure 1: Robustness and reproducibility of skin metatranscriptomes across 1206 different skin sites . (A) Boxplot showing the fraction of non -ribosomal RNA (non 1207 rRNA) reads with and without experimental rRNA depletion during library preparation 1208 for the pilot cohort . (B) Boxplot of Sorensen Similarity Ind ices (1 – Bray Curtis 1209 dissimilarity) computed from species -level relative abundances for the pilot cohort . 1210 Pairwise similarities were computed between samples in three different categories. 1211 “Technical replicates ” refer to different RNA -seq librar ies prepared and sequenc ed 1212 from the same samples. “Within individuals” refer to data from re -sampling the same 1213 individuals across three consecutive days, while “between individuals” represents 1214 inter-individual variability in our cohort. (C) Same as B but showing Pearson correlation 1215 of gene expression signatures in the pilot cohort . (D) Schematic of the full cohort 1216 comprising 27 healthy adult volunteers and 5 different skin sites and the data analysis 1217 workflow. Boxplot showing the number of metatranscriptomic reads before 1218 computational removal of host and rRNA reads (red), number of non -human, non -1219 rRNA reads after computational filtering (green) and number of microbial reads 1220 (comprising bacteria, viruses, fungi and archaea) after de -duplication (blue). (E) 1221 Boxplot showing the fraction of non -ribosomal RNA (non rRNA) reads with 1222 experimental rRNA depletion during library preparation for the full cohort. (F) Coverage 1223 plots showing the distribution of reads over bacterial gene bodies in samples belonging 1224 to the full cohort . Reads were mapped to bacterial pangenomes and samples with 1225 either >50,000 or ≤50,000 mapped bacteria l reads are coloured differently. (G) 1226 Coverage plots showing the distribution of reads over fungal gene bodies for libraries 1227 (full cohort) with ≥500,000 fungal reads. (H) Percentage of reads mapped to different 1228 regions (intergenic, intronic and exons) of the Malassezia globosa genome. “This 1229 study” refers to samples from the full cohort. “Wu et al 2015” refers to an external RNA-1230 seq dataset from cultured Malassezia globosa isolates27. 1231 Figure 2: Niche specific signatures of core active species and transcribed 1232 functions. (A) Bar plots showing median relative abundances (total sum scaled) of 1233 reads assigned to various skin commensals across different skin sites presented in 1234 log2 scale. Abundances for metagenomes and metatranscriptomes are shown. (B) 1235 Bubble plot showing core, common and variable components of skin metagenomes 1236 (DNA) and metatranscriptomes (RNA) across different body sites . A species was 1237 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint called present in a metagenome if it had ≥0.1% relative abundance. A species was 1238 called present in a metatranscriptome if it had ≥0.1% relative abundance and was also 1239 detected in the metagenome. Core species in metagenomes and metatranscriptomes 1240 were defined as those present across >75% of samples at a given skin site. Common 1241 species for a skin site were defined as those present between 50% and 75% of 1242 samples. Variable species for a given skin site refer to other species that do not fall 1243 into the previous two categories, but which are present in ≥3 individuals. For RNA, 1244 bubbles are also shaded according to median transcriptional activity , defined as the 1245 ratio of normalized RNA counts of a species to that of their DNA. For visual clarity, 1246 only species which were core in at least one skin site are represented here. (C) 1247 Scatterplot of mean beta diversity (1 - Bray-Curtis dissimilarity) against mean alpha 1248 diversity (Simpson index) of core microbial pathways. Core microbial pathway 1249 expression was computed using HUMAnN3. Core microbial pathways were defined 1250 as those which were present (non-zero expression) at a skin site in >75% of individuals 1251 and with <25% unclassified reads at species -level. (D) Stacked bar plots for species 1252 level pathway contributions at RNA level, estimated with HUMAnN3 for Staphyloferrin 1253 A biosynthesis, and with Kraken2 for c ommunity level relative abundances at RNA 1254 level. For all sub-figures, Sc, Ch, Ac, Vf and Tw indicate the skin sites scalp, cheek, 1255 antecubital fossa, volar forearm and toe web respectively. 1256 Figure 3: Species-level differential enrichment of metabolic pathways in various 1257 in vivo and in vitro growth conditions. (A) Volcano plot of differentially expressed 1258 genes for Malassezia restricta colonizing scalp versus cheek samples (≥200,000 M. 1259 restricta reads per library) . Genes involved in peroxisomal activity or ether 1260 lipid/glycerophospholipid/sphingolipid metabolism are coloured blue or yellow 1261 respectively. The genes PLCN_1-7 encode secreted phosopholipase C enzymes. The 1262 genes POX2, POT1, PEX11B and ACOT8 encode proteins for peroxisomal function. 1263 (B) Bar plots of Normalized Enrichment Scores (NES) from Gene Set Enrichment 1264 Analysis for Malassezia restricta genes differentially expressed in scalp versus cheek 1265 samples. (C) Bar plots showing the number of differentially expressed genes between 1266 scalp and cheek sites for different species. (D) Schematic showing experimental setup 1267 for in vitro cultures of Staphylococcus epidermidis under different growth conditions 1268 and stress exposures. (E) Left: Principal component analysis (PCA) plot for S. 1269 epidermidis gene expression profiles from various in vitro and in vivo samples after 1270 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint batch correction using Limma. Right: Bar plots showing the number of differentially 1271 expressed genes for S. epidermidis between in vivo transcriptomes and 1272 transcriptomes for different in vitro conditions. Sebaceous sites refer to both scalp and 1273 cheek samples with ≥200,000 S. epidermidis reads per library (n=6). Toe webs refer 1274 to toe web samples with ≥200,000 S. epidermidis reads per library (n=12). 1275 Figure 4: Transcriptome-aware, genome scale metabolic models reveal 1276 metabolic dependencies in various in vivo and in vitro growth conditions. (A) 1277 Bar plots showing mean flux values of various reactions associated with propionate 1278 production during cysteine and methionine metabolism in different in vivo and in vitro 1279 growth conditions for Staphylococcus epidermidis. Flux values were estimated using 1280 transcriptome-aware genome scale metabolic models for S. epidermidis. EX_ppae(e) 1281 represents the reaction exporting propionate out of bacterial cells. All reaction fluxes 1282 are significantly higher for in vivo relative to in vitro conditions (GLM adjusted p -1283 value<0.001). (B) Same as A, except showing various reactions involved in pyruvate 1284 generation during glycolysis. All reaction fluxes are significantly different between in 1285 vivo and in vitro conditions. Reaction names highlighted in red are alternative means 1286 of generating pyruvate that have higher mean fluxes for in vivo versus in vitro 1287 conditions, consistent with in vivo metabolic dependency. (C) Bar plots showing mean 1288 flux values of various reactions associated with propionate production during cysteine 1289 and methionine metabolism in different in vivo growth conditions (Scalp or Cheek) for 1290 Cutibacterium acnes. Flux values were estimated using transcriptome-aware genome 1291 scale metabolic models for C. acnes. (D) Same as C, except showing various 1292 reactions involved in histidine and proline/glutamate metabolism. EX_his_L(e) and 1293 EX_pro_L(e) refer to transport reactions exporting histidine and proline respectively. 1294 Figure 5: Inferring host-microbe and microbe-microbe interactions from 1295 metatranscriptomes. (A) Heatmap of various classes of antimicrobial genes with 1296 detected expression (≥5 reads) in individual metatranscriptomes grouped by skin site. 1297 Blue tiles denote detected expression in a given metatranscriptome. Red labels denote 1298 genes with site -specific enrichments based on Fisher’s exact test. Bacterial genes 1299 were grouped into categories of antimicrobial features based on information from 1300 profile Hidden Markov Models downloaded from NCBI. (B) Heatmap of various classes 1301 of antimicrobial genes with expression associated with species -level classified 1302 metatranscriptomic reads. The number of metatranscriptomes with detected 1303 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint expression (≥5 reads) is shown in the figure . (C) Overlaid protein structures of the 1304 Malassezia restricta protein DNF11_2196 and a Streptomyces papain inhibitor. 1305 Structures were obtained from the publicly available Alpha Fold database. Structural 1306 similarity scores from Foldseek (TM-score and Root Mean Square Deviation) and Dali 1307 (Z-score) are shown. Primary amino acid sequence identity between the two proteins 1308 is also given. (D) Left: Boxplots of pro -IL-1B levels from human keratinocytes as 1309 measured by HiBiT and averaged across multiple strains of Staphylococcus capitis 1310 (n=4), Staphylococcus epidermidis (n=4), Staphylococcus hominis (n=8) and 1311 Staphylococcus warneri (n=4). Data are from two biological repeats, each comprising 1312 of three technical replicates. A Kruskal -Wallis test was conducted to assess whether 1313 levels differed significantly between groups (Kruskal-Wallis chi-squared=19.854, df=5, 1314 p-value<0.0014) and the adjusted p -values from post -hoc Dunn tests are shown. 1315 Right: Boxplots of cleaved IL -1B levels from human keratinocytes as measured by 1316 ELISA and averaged across multiple strains of Staphylococcus capitis , 1317 Staphylococcus epidermidis , Staphylococcus hominis and Staphylococcus warneri . 1318 Data are from two biological repeats, each comprising of three technical replicates. A 1319 Kruskal-Wallis test was conducted to assess whether levels differed significantly 1320 between groups (Kruskal-Wallis chi-squared=33.864, df=5, p-value<2.6e-06) and the 1321 adjusted p-values from post-hoc Dunn tests are shown. 1322 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint

1323 1. Byrd, A. L., Belkaid, Y. & Segre, J. A. The human skin microbiome. Nature Reviews 1324 Microbiology vol. 16 143–155 Preprint at https://doi.org/10.1038/nrmicro.2017.157 1325 (2018). 1326 2. Claesen, J. et al. A cutibacterium acnes antibiotic modulates human skin microbiota 1327 composition in hair follicles. Sci Transl Med 12, (2020). 1328 3. Tay, A. S. L. et al. Atopic dermatitis microbiomes stratify into ecologic dermotypes 1329 enabling microbial virulence and disease severity. Journal of Allergy and Clinical 1330 Immunology (2020) doi:10.1016/j.jaci.2020.09.031. 1331 4. Oh, J. et al. Biogeography and individuality shape function in the human skin 1332 metagenome. Nature 514, 59–64 (2014). 1333 5. Uberoi, A. et al. Commensal microbiota regulates skin barrier function and repair via 1334 signaling through the aryl hydrocarbon receptor. Cell Host Microbe 29, 1235-1248.e8 1335 (2021). 1336 6. Deng, L. et al. S. aureus drives itch and scratch-induced skin damage through a V8 1337 protease-PAR1 axis. Cell 186, 5375-5393.e25 (2023). 1338 7. Li, H. et al. Skin Commensal Malassezia globosa Secreted Protease Attenuates 1339 Staphylococcus aureus Biofilm Formation. Journal of Investigative Dermatology 138, 1340 1137–1145 (2018). 1341 8. Williams, M. R. et al. Quorum sensing between bacterial species on the skin protects 1342 against epidermal injury in atopic dermatitis. Sci Transl Med 11, 8329 (2019). 1343 9. Franzosa, E. A. et al. Species-level functional profiling of metagenomes and 1344 metatranscriptomes. Nat Methods 15, 962–968 (2018). 1345 10. Poh, S. E. et al. Expression of Staphylococcus aureus Virulence Factors in Atopic 1346 Dermatitis. JID Innov 2, 100130 (2022). 1347 11. Chua, W. et al. A Malassezia pseudoprotease dominates the secreted hydrolase 1348 landscape and is a potential allergen on skin. Biochimie 216, 181–193 (2024). 1349 12. Chaves-Moreno, D. et al. Exploring the transcriptome of Staphylococcus aureus in its 1350 natural niche. Scientific Reports 2016 6:1 6, 1–11 (2016). 1351 13. Abu-Ali, G. S. et al. Metatranscriptome of human faecal microbial communities in a 1352 cohort of adult men. Nature Microbiology 2018 3:3 3, 356–366 (2018). 1353 14. Salazar, G. et al. Gene Expression Changes and Community Turnover Differentially 1354 Shape the Global Ocean Metatranscriptome. Cell 179, 1068-1083.e21 (2019). 1355 15. Whitman, W. B., Coleman, D. C. & Wiebe, W. J. Prokaryotes: The unseen majority. 1356 Proc Natl Acad Sci U S A 95, 6578 (1998). 1357 16. Kang, D., Shi, B., Erfe, M. C., Craft, N. & Li, H. Vitamin B12 modulates the 1358 transcriptome of the skin microbiota in acne pathogenesis. Sci Transl Med 7, 1359 293ra103 (2015). 1360 17. Kong, H. H. et al. Performing Skin Microbiome Research: A Method to the Madness. J 1361 Invest Dermatol 137, 561–568 (2017). 1362 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint 18. Santiago-Rodriguez, T. M., Le François, B., Macklaim, J. M., Doukhanine, E. & 1363 Hollister, E. B. The Skin Microbiome: Current Techniques, Challenges, and Future 1364 Directions. Microorganisms 11, (2023). 1365 19. Li, Z. et al. Characterization of the human skin resistome and identification of two 1366 microbiota cutotypes. Microbiome 9, 1–18 (2021). 1367 20. Beghini, F. et al. Integrating taxonomic, functional, and strain-level profiling of diverse 1368 microbial communities with bioBakery 3. Elife 10, (2021). 1369 21. Tan, C. C. S. et al. No evidence for a common blood microbiome based on a 1370 population study of 9,770 healthy humans. Nat Microbiol 8, (2023). 1371 22. de Goffau, M. C. et al. Human placenta has no microbiome but can contain potential 1372 pathogens. Nature 572, 329–334 (2019). 1373 23. Breitwieser, F. P., Baker, D. N. & Salzberg, S. L. KrakenUniq: Confident and fast 1374 metagenomics classification using unique k-mer counts. Genome Biol 19, 1–10 1375 (2018). 1376 24. Wood, D. E., Lu, J. & Langmead, B. Improved metagenomic analysis with Kraken 2. 1377 Genome Biol 20, (2019). 1378 25. Franzosa, E. A. et al. Relating the metatranscriptome and metagenome of the human 1379 gut. Proc Natl Acad Sci U S A 111, E2329–E2338 (2014). 1380 26. Deng, Z. L., Münch, P . C., Mreches, R. & McHardy, A. C. Rapid and accurate 1381 identification of ribosomal RNA sequences via deep learning. Nucleic Acids Res 50, 1382 e60 (2022). 1383 27. Wu, G. et al. Genus-Wide Comparative Genomics of Malassezia Delineates Its 1384 Phylogeny, Physiology, and Niche Adaptation on Human Skin. PLoS Genet 11, 1385 1005614 (2015). 1386 28. Vijaya Chandra, S. H., Srinivas, R., Dawson, T. L. & Common, J. E. Cutaneous 1387 Malassezia: Commensal, Pathogen, or Protector? Front Cell Infect Microbiol 10, 1388 614446 (2021). 1389 29. Chong, C. E., Bengtsson, R. J. & Horsburgh, M. J. Comparative genomics of 1390 Staphylococcus capitis reveals species determinants. Front Microbiol 13, 1005949 1391 (2022). 1392 30. Swaney, M. H., Nelsen, A., Sandstrom, S. & Kalan, L. R. Sweat and Sebum 1393 Preferences of the Human Skin Microbiota. Microbiol Spectr 11, (2023). 1394 31. Prohic, A., Jovovic Sadikovic, T., Krupalija-Fazlic, M. & Kuskunovic-Vlahovljak, S. 1395 Malassezia species in healthy skin and in dermatological conditions. Int J Dermatol 1396 55, 494–504 (2016). 1397 32. Magiatis, P . et al. Malassezia yeasts produce a collection of exceptionally potent 1398 activators of the ah (dioxin) receptor detected in diseased human skin. Journal of 1399 Investigative Dermatology 133, 2023–2030 (2013). 1400 33. Murphy, G. R. et al. Relationships between electrolyte and amino acid compositions in 1401 sweat during exercise suggest a role for amino acids and K+ in reabsorption of Na+ 1402 and Cl- from sweat. PLoS One 14, (2019). 1403 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint 34. Imi, Y . & Shibata, K. Nutritional Factors That Affect the Formation of 5-Aminolevulinic 1404 Acid, a Key Intermediate of Heme Biosynthesis. J Nutr Sci Vitaminol (Tokyo) 67, 339–1405 350 (2021). 1406 35. Leasure, C. S., Grunenwald, C. M., Choby, J. E., Sauer, J. D. & Skaar, E. P. 1407 Maintenance of heme homeostasis in Staphylococcus aureus through post-1408 translational regulation of glutamyl-tRNA reductase. J Bacteriol 205, (2023). 1409 36. Thiele, J. J., Weber, S. U. & Packer, L. Sebaceous Gland Secretion is a Major 1410 Physiologic Route of Vitamin E Delivery to Skin. Journal of Investigative Dermatology 1411 113, 1006–1010 (1999). 1412 37. Leung, M. H. Y . et al. Skin microbiome differentiates into distinct cutotypes with unique 1413 metabolic functions upon exposure to polycyclic aromatic hydrocarbons. Microbiome 1414 11, 1–14 (2023). 1415 38. Zheng, Y . et al. Commensal Staphylococcus epidermidis contributes to skin barrier 1416 homeostasis by generating protective ceramides. Cell Host Microbe 30, 301-313.e9 1417 (2022). 1418 39. Jourdain, R. et al. Exploration of scalp surface lipids reveals squalene peroxide as a 1419 potential actor in dandruff condition. Arch Dermatol Res 308, 153 (2016). 1420 40. Baker, L. B. Physiology of sweat gland function: The roles of sweating and sweat 1421 composition in human health. Temperature: Multidisciplinary Biomedical Journal 6, 1422 211 (2019). 1423 41. Wang, X. et al. PhoU2 but not PhoU1 as an important regulator of biofilm formation 1424 and tolerance to multiple stresses by participating in various fundamental metabolic 1425 processes in Staphylococcus epidermidis. J Bacteriol 199, (2017). 1426 42. Avican, K. et al. RNA atlas of human bacterial pathogens uncovers stress dynamics 1427 linked to infection. Nature Communications 2021 12:1 12, 1–14 (2021). 1428 43. Palmgren, M. P-type ATPases: Many more enigmas left to solve. Journal of Biological 1429 Chemistry 299, 105352 (2023). 1430 44. Cohn, J. R. & Emmett, E. A. The excretion of trace metals in human sweat. Ann Clin 1431 Lab Sci 8, 270–5 (1978). 1432 45. Kuan, W. H., Chen, Y . L. & Liu, C. L. Excretion of Ni, Pb, Cu, As, and Hg in Sweat 1433 under Two Sweating Conditions. International Journal of Environmental Research and 1434 Public Health 2022, Vol. 19, Page 4323 19, 4323 (2022). 1435 46. Chandrangsu, P ., Rensing, C. & Helmann, J. D. Metal homeostasis and resistance in 1436 bacteria. Nature Reviews Microbiology 2017 15:6 15, 338–350 (2017). 1437 47. Jenior, M. L., Moutinho, T. J., Dougherty, B. V. & Papin, J. A. Transcriptome-guided 1438 parsimonious flux analysis improves predictions with metabolic networks in complex 1439 environments. PLoS Comput Biol 16, e1007099 (2020). 1440 48. Xiao, X. et al. The role of short-chain fatty acids in inflammatory skin diseases. Front 1441 Microbiol 13, 1083432 (2023). 1442 49. Gonzalez-Garcia, R. A., McCubbin, T., Turner, M. S., Nielsen, L. K. & Marcellin, E. 1443 Engineering Escherichia coli for propionic acid production through the Wood–1444 Werkman cycle. Biotechnol Bioeng 117, 167–183 (2020). 1445 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint 50. Walker, M. C. & van der Donk, W. A. The many roles of glutamate in metabolism. J 1446 Ind Microbiol Biotechnol 43, 419–430 (2016). 1447 51. Heilbronner, S., Krismer, B., Brötz-Oesterhelt, H. & Peschel, A. The microbiome-1448 shaping roles of bacteriocins. Nature Reviews Microbiology 2021 19:11 19, 726–739 1449 (2021). 1450 52. Nakazono, K. et al. Complete sequences of epidermin and nukacin encoding 1451 plasmids from oral-derived Staphylococcus epidermidis and their antibacterial activity. 1452 PLoS One 17, (2022). 1453 53. Ebner, P . et al. Lantibiotic production is a burden for the producing staphylococci. 1454 Scientific Reports 2018 8:1 8, 1–11 (2018). 1455 54. Atanasova, N. S., Pietilä, M. K. & Oksanen, H. M. Diverse antimicrobial interactions of 1456 halophilic archaea and bacteria extend over geographical distances and cross the 1457 domain barrier. Microbiologyopen 2, 811–825 (2013). 1458 55. Makarova, K. S. et al. Antimicrobial peptides, polymorphic toxins, and self-nonself 1459 recognition systems in archaea: An untapped armory for intermicrobial conflicts. mBio 1460 10, (2019). 1461 56. Ziegenhain, C. et al. Comparative Analysis of Single-Cell RNA Sequencing Methods. 1462 Mol Cell 65, 631-643.e4 (2017). 1463 57. Wu, A. R. et al. Quantitative assessment of single-cell RNA-sequencing methods. 1464 Nature Methods 2013 11:1 11, 41–46 (2013). 1465 58. Sølberg, J. et al. The stratum corneum transcriptome in atopic dermatitis can be 1466 assessed by tape stripping. J Dermatol Sci 101, 14–21 (2021). 1467 59. Amorim, C. F. et al. Variable gene expression and parasite load predict treatment 1468 outcome in cutaneous leishmaniasis. Sci Transl Med 11, (2019). 1469 60. Hänzelmann, S., Castelo, R. & Guinney, J. GSVA: Gene set variation analysis for 1470 microarray and RNA-Seq data. BMC Bioinformatics 14, 1–15 (2013). 1471 61. O’Shea, J. J. et al. Signal transduction and Th17 cell differentiation. Microbes Infect 1472 11, 599–611 (2009). 1473 62. Tham, K. C. et al. Distinct skin microbiome community structures in congenital 1474 ichthyosis. Br J Dermatol 187, 557 (2022). 1475 63. Jumper, J. et al. Highly accurate protein structure prediction with AlphaFold. Nature 1476 2021 596:7873 596, 583–589 (2021). 1477 64. van Kempen, M. et al. Fast and accurate protein structure search with Foldseek. 1478 Nature Biotechnology 2023 42:2 42, 243–246 (2023). 1479 65. Zindel, S. et al. The Papain Inhibitor (SPI) of Streptomyces mobaraensis Inhibits 1480 Bacterial Cysteine Proteases and Is an Antagonist of Bacterial Growth. Antimicrob 1481 Agents Chemother 57, 3388 (2013). 1482 66. Smith, T., Heger, A. & Sudbery, I. UMI-tools: Modeling sequencing errors in Unique 1483 Molecular Identifiers to improve quantification accuracy. Genome Res 27, 491–499 1484 (2017). 1485 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint 67. Hall, J. B. et al. Isolation and Identification of the Follicular Microbiome: Implications 1486 for Acne Research. Journal of Investigative Dermatology 138, 2033–2040 (2018). 1487 68. Nestle, F. O., Di Meglio, P., Qin, J. Z. & Nickoloff, B. J. Skin immune sentinels in 1488 health and disease. Nature Reviews Immunology 2009 9:10 9, 679–691 (2009). 1489 69. Conwill, A. et al. Anatomy promotes neutral coexistence of strains in the human skin 1490 microbiome. Cell Host Microbe 30, 171-182.e7 (2022). 1491 70. Magiatis, P . et al. Malassezia yeasts produce a collection of exceptionally potent 1492 activators of the Ah (dioxin) receptor detected in diseased human skin. J Invest 1493 Dermatol 133, 2023–2030 (2013). 1494 71. Plotkin, L. I., Mathov, I., Squiquera, L. & Leoni, J. Arachidonic acid released from 1495 epithelial cells by Malassezia furfur phospholipase A2: a potential pathophysiologic 1496 mechanism. Mycologia 90, 163–169 (1998). 1497 72. Potter, A. D. et al. Host nutrient milieu drives an essential role for aspartate 1498 biosynthesis during invasive Staphylococcus aureus infection. Proc Natl Acad Sci U S 1499 A 117, 12394–12401 (2020). 1500 73. DeMars, Z., Singh, V. K. & Bose, J. L. Exogenous fatty acids remodel staphylococcus 1501 aureus lipid composition through fatty acid kinase. J Bacteriol 202, (2020). 1502 74. Onyango, L. A. & Alreshidi, M. M. Adaptive Metabolism in Staphylococci: Survival and 1503 Persistence in Environmental and Clinical Settings. J Pathog 2018, 1092632 (2018). 1504 75. Perez, K. & Patel, R. Survival of Staphylococcus epidermidis in Fibroblasts and 1505 Osteoblasts. Infect Immun 86, (2018). 1506 76. Naik, S. et al. Commensal–dendritic-cell interaction specifies a unique protective skin 1507 immune signature. Nature 2015 520:7545 520, 104–108 (2015). 1508 77. Lam, T. H. et al. Understanding the microbial basis of body odor in pre-pubescent 1509 children and teenagers. Microbiome 6, 1–14 (2018). 1510 78. Jenior, M. L., Moutinho, T. J., Dougherty, B. V. & Papin, J. A. Transcriptome-guided 1511 parsimonious flux analysis improves predictions with metabolic networks in complex 1512 environments. PLoS Comput Biol 16, e1007099 (2020). 1513 79. Bordbar, A., Monk, J. M., King, Z. A. & Palsson, B. O. Constraint-based models 1514 predict metabolic and associated cellular functions. Nature Reviews Genetics 2014 1515 15:2 15, 107–120 (2014). 1516 80. Kim, J., Kim, G. L., Norambuena, J., Boyd, J. M. & Parker, D. Impact of the pentose 1517 phosphate pathway on metabolism and pathogenesis of Staphylococcus aureus. 1518 PLoS Pathog 19, e1011531 (2023). 1519 81. Sanford, J. A., O’Neill, A. M., Zouboulis, C. C. & Gallo, R. L. Short-Chain Fatty Acids 1520 from Cutibacterium acnes Activate Both a Canonical and Epigenetic Inflammatory 1521 Response in Human Sebocytes. The Journal of Immunology 202, 1767–1776 (2019). 1522 82. Sashihara, T. et al. A Novel Lantibiotic, Nukacin ISK-1, of Staphylococcus warneri ISK-1523 1: Cloning of the Structural Gene and Identification of the Structure. Biosci Biotechnol 1524 Biochem 64, 2420–2428 (2000). 1525 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint 83. DI Tommaso, P. et al. Nextflow enables reproducible computational workflows. Nature 1526 Biotechnology 2017 35:4 35, 316–319 (2017). 1527 84. Chen, S., Zhou, Y., Chen, Y . & Gu, J. fastp: an ultra-fast all-in-one FASTQ 1528 preprocessor. Bioinformatics 34, i884–i890 (2018). 1529 85. Li, H. & Durbin, R. Fast and accurate short read alignment with Burrows–Wheeler 1530 transform. Bioinformatics 25, 1754 (2009). 1531 86. Li, H. et al. The Sequence Alignment/Map format and SAMtools. Bioinformatics 25, 1532 2078–2079 (2009). 1533 87. Dobin, A. et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29, 15 1534 (2013). 1535 88. Lu, J., Breitwieser, F. P ., Thielen, P . & Salzberg, S. L. Bracken: Estimating species 1536 abundance in metagenomics data. PeerJ Comput Sci 2017, e104 (2017). 1537 89. Janssens, Y . et al. Disbiome database: Linking the microbiome to disease. BMC 1538 Microbiol 18, 1–6 (2018). 1539 90. Yao, G. et al. MicroPhenoDB Associates Metagenomic Data with Pathogenic 1540 Microbes, Microbial Core Genes, and Human Disease Phenotypes. Genomics 1541 Proteomics Bioinformatics 18, 760–772 (2020). 1542 91. Friedman, J. & Alm, E. J. Inferring Correlation Networks from Genomic Survey Data. 1543 PLoS Comput Biol 8, e1002687 (2012). 1544 92. Tan, C. C. S. et al. No evidence for a common blood microbiome based on a 1545 population study of 9,770 healthy humans. Nature Microbiology 2023 8:5 8, 973–985 1546 (2023). 1547 93. Langmead, B. & Salzberg, S. L. Fast gapped-read alignment with Bowtie 2. Nat 1548

9, 357 (2012). 1549 94. Suzek, B. E., Wang, Y., Huang, H., McGarvey, P . B. & Wu, C. H. UniRef clusters: a 1550 comprehensive and scalable alternative for improving sequence similarity searches. 1551 Bioinformatics 31, 926 (2015). 1552 95. Buchfink, B., Reuter, K. & Drost, H. G. Sensitive protein alignments at tree-of-life 1553 scale using DIAMOND. Nature Methods 2021 18:4 18, 366–368 (2021). 1554 96. Cantalapiedra, C. P., Hern ̗ andez-Plaza, A., Letunic, I., Bork, P. & Huerta-Cepas, J. 1555 eggNOG-mapper v2: Functional Annotation, Orthology Assignments, and Domain 1556 Prediction at the Metagenomic Scale. Mol Biol Evol 38, 5825–5829 (2021). 1557 97. Huerta-Cepas, J. et al. eggNOG 5.0: a hierarchical, functionally and phylogenetically 1558 annotated orthology resource based on 5090 organisms and 2502 viruses. Nucleic 1559 Acids Res 47, D309–D314 (2019). 1560 98. Consortium, T. G. O. et al. The Gene Ontology knowledgebase in 2023. Genetics 224, 1561 (2023). 1562 99. Ashburner, M. et al. Gene Ontology: tool for the unification of biology. Nat Genet 25, 1563 25 (2000). 1564 100. Kanehisa, M. & Goto, S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic 1565 Acids Res 28, 27–30 (2000). 1566 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint 101. Zhang, Y . et al. ProPan: a comprehensive database for profiling prokaryotic pan-1567 genome dynamics. Nucleic Acids Res 51, D767–D776 (2023). 1568 102. Patro, R., Duggal, G., Love, M. I., Irizarry, R. A. & Kingsford, C. Salmon: fast and bias-1569 aware quantification of transcript expression using dual-phase inference. Nat Methods 1570 14, 417 (2017). 1571 103. Love, M. I., Huber, W. & Anders, S. Moderated estimation of fold change and 1572 dispersion for RNA-seq data with DESeq2. Genome Biol 15, 1–21 (2014). 1573 104. Ritchie, M. E. et al. limma powers differential expression analyses for RNA-1574 sequencing and microarray studies. Nucleic Acids Res 43, e47–e47 (2015). 1575 105. Yu, G., Wang, L. G., Han, Y . & He, Q. Y. clusterProfiler: an R Package for Comparing 1576 Biological Themes Among Gene Clusters. OMICS 16, 284 (2012). 1577 106. Eddy, S. R. Accelerated Profile HMM Searches. PLoS Comput Biol 7, (2011). 1578 107. Edgar, R. C. Muscle5: High-accuracy alignment ensembles enable unbiased 1579 assessments of sequence homology and phylogeny. Nat Commun 13, (2022). 1580 108. Teufel, F. et al. SignalP 6.0 predicts all five types of signal peptides using protein 1581 language models. Nature Biotechnology 2022 40:7 40, 1023–1025 (2022). 1582 109. Fernandes, A. D. et al. Unifying the analysis of high-throughput sequencing datasets: 1583 characterizing RNA-seq, 16S rRNA gene sequencing and selective growth 1584 experiments by compositional data analysis. Microbiome 2, 15 (2014). 1585 110. Magnúsdóttir, S. et al. Generation of genome-scale metabolic reconstructions for 773 1586 members of the human gut microbiota. Nature Biotechnology 2016 35:1 35, 81–89 1587 (2016). 1588 111. Ebrahim, A., Lerman, J. A., Palsson, B. O. & Hyduke, D. R. COBRApy: COnstraints-1589 Based Reconstruction and Analysis for Python. BMC Syst Biol 7, 1–6 (2013). 1590 112. Jenior, M. L., Moutinho, T. J., Dougherty, B. V. & Papin, J. A. Transcriptome-guided 1591 parsimonious flux analysis improves predictions with metabolic networks in complex 1592 environments. PLoS Comput Biol 16, e1007099 (2020). 1593 1594 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint Species leveltechnical replicates within individuals between individuals Sorensen Similarity Index 0.9 0.6 0.3 0.0 Gene family level A C DB Ac Sc Ch Sc Ch Ac Vf Tw 1.0 0.5 1.0 0.5 5' 3' Vf Tw Map to bacterial pangenomes Map to M. globosa/M. restricta CDS 5' 3'CDS normalized coverage normalized coverage Pearson correlation <=50k mapped reads 1.0 0.5 1.0 0.5 1.0 0.5 (n=23) (n=19) (n=16) (n=7)(n=18) (n=18) (n=24) (n=18) (n=18) (n=1) 1.0 0.5 1.0 0.6 0.9 0.3 0.9 0.3 1.0 0.5 % mapping to genomic regions of M. globosa Intergenic Intronic mRNA 80 60 40 20 0 percentage mapped thisstudy Wu etal 2015 thisstudy Wu etal 2015 thisstudy Wu etal 2015 QC n = 27 b) Cheek (Ch) e) T oe webs (Tw) Sebaceous Moist Dry Full cohort d) Volar forearm (Vf) a) Scalp (Sc) c) Antecubital fossa (Ac) Annotate Kitome removal T axonomic classificationSkin microbial reads Functional inference n=19 Ac n=18 Vf n=18 T w n=23n=24 Sc Ch Sc Ch Ac Vf T w 0.00 0.25 0.50 0.75 1.00 p < 2.2e−06 p < 6e−12 technical replicates within individuals between individuals p < 1.8e−12 p < 6.8e−09 No rRNA depletion Yes p = 0.014 Figure 1 75 100 50 25 80 100 60 40 % non rRN A reads% non rRN A reads All reads non-human, non-rRNAs microbial reads de-duplicated 6.0 6.5 7.0 7.5 8.0log10 read pairs Metatranscriptomes (n=102) E F G H >50k mapped reads .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint Sc Ch Ac Vf T w Sc Ch Ac Vf T w core (>75%) common (>50%) variable (≥3 individuals) Transcriptional activity log2 median (RNA/DNA) 5 0 10 -5 Prevalence at site (%) 100 80 60 40 20 Corynebacterium kefirresidentii Corynebacterium segmentosum Corynebacterium tuberculostearicum Corynebacterium ureicelerivorans Cutibacterium acnes Cutibacterium granulosum Cutibacterium modestum Deinococcus wulumuqiensis Malassezia arunalokei Malassezia dermatis Malassezia furfur Malassezia globosa Malassezia restricta Malassezia sympodialis Moraxella osloensis Rothia mucilaginosa Staphylococcus capitis Staphylococcus caprae Staphylococcus epidermidis Staphylococcus haemolyticus Staphylococcus hominis Staphylococcus warneri Figure 2 0.1 0.2 0.4 0.6 0.8 0.2 0.3 0.4 0.5 0.1 0.2 0.3 0.4 0.5 0.1 0.2 0.3 0.4 0.5 0.1 0.2 0.3 0.4 0.5 0.1 0.2 0.3 0.4 0.5 Mean alpha diversity Mean beta diversity Sc Ch Ac Vf Tw High alpha, High beta Low alpha, High beta Low alpha, Low beta Community taxonomic diversity RNADNA Scalp log2 (median rel. abundance) −3 0 3 6 −2.5 0.0 2.5 5.0 −2.5 0.0 2.5 5.0 Cheek Antecubital fossa Metagenomes Metatranscriptomes Volar forearm T oe web Staphylococcus capitis Staphylococcus epidermidis Staphylococcus hominis Malassezia globosa Malassezia restricta Cutibacterium acnes Staphylococcus capitis Staphylococcus epidermidis Staphylococcus hominis Malassezia globosa Malassezia restricta Cutibacterium acnes Staphylococcus capitis Staphylococcus epidermidis Staphylococcus hominis Malassezia globosa Malassezia restricta Cutibacterium acnes Staphylococcus capitis Staphylococcus epidermidis Staphylococcus hominis Malassezia globosa Malassezia restricta Cutibacterium acnes Staphylococcus capitis Staphylococcus epidermidis Staphylococcus hominis Malassezia globosa Malassezia restricta Cutibacterium acnes −2 0 2 4 −2.5 0.0 2.5 Staphylococcus capitis Staphylococcus caprae Staphylococcus epidermidis Staphylococcus haemolyticus Staphylococcus hominis Staphylococcus warneri others unclassified at species level L-ornithine Staphyloferrin A 0 25 50 75 100 0 25 50 75 100 Pathway (M00876) T oe web communitySMT001_Tw SMT002_Tw SMT003_Tw SMT008_Tw SMT012_Tw SMT014_Tw SMT015_Tw SMT016_Tw SMT017_Tw SMT018_Tw SMT019_Tw SMT020_Tw SMT021_Tw SMT023_Tw SMT026_Tw SMT001_Tw SMT002_Tw SMT003_Tw SMT008_Tw SMT012_Tw SMT014_Tw SMT015_Tw SMT016_Tw SMT017_Tw SMT018_Tw SMT019_Tw SMT020_Tw SMT021_Tw SMT023_Tw SMT026_Tw Low alpha diversity High beta diversity Pathway relativecontribution (RNA) Species relativeabundance (RNA) B C D A .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint Figure 3 −2 −1 0 31 2 Normalized Enrichment Score STT4 PEX11B ACOT8 POX2 PLCN_5 PLCN_6 PLCN_3PLCN_4 LACS7 POT1 PLCN_1 PLCN_2 PLCN_7 PIK1 Scalp Scalp Oily Cheek Cheek Oily Inositol phosphate metabolism Quorum sensing Glycerophospholipid metabolism Ether lipid metabolism Cell cycle − yeast Microbial metabolism diverse environments Protein processing in ER Fructose & mannose metabolism Proteasome Peroxisome Ribosome Peroxisome KEGG Pathways -Log10 (p-value) Log2 fold change Ether lipid/Glycerophospholipid/Sphingolipid metabolism Malassezia restricta Scalp vs CheekMalassezia restricta Differentially expressed genes 300 200 100 0 C. acnesM. globosaM. restricta upregulated in scalp downregulated in scalp −50 −25 0 25 50 −50 −25 0 25 50 PC1: 66% variance PC2: 8% variance Log phase Osmotic stress Stationary phase in vitro Sebaceous site T oe Web in vivo upregulated in vivo downregulated in vivo 0 250 500 750 1000 log phase osmotic stress stationary phase log phase osmotic stress stationary phase Sebaceous site vs. T oe web vs. S. epidermidis in vivo versus in vitro Differentially expressed genes 1) Log phase in vitro dataset other published studies 2) Osmotic stress 3) Stationary phase sub- culture S. epidermidis overnight cultures log phase cultures 0.5M NaCl 16h 30 -5.0 -2.5 2.50.0 20 10 0 D E A B C .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint A B C D Figure 4 Cys & Met metabolism 2-Oxobutanoate L-Cystathionine Propanoyl phosphate Propionate OBTFL & PTA2 PPAKr CYSTGL Propionate metabolism Staphylococcus epidermidis Staphylococcus epidermidis Glycolytic pathway Glucose F6P G3P DHAP S7P Pyruvate FDP Phosphotransferase system Various enzymes L-Lactate dehydrogenase (L_LACD2) Transaldolase (TALA) PFK (PFK_3) S1,7P G3P lyase (FBA3) L-Lactate Glucose-6-phosphate isomerase Pentose Phosphate Pathway 0 2 4 6 8 CYSTGL OBTFL PPAKr PTA2 EX_ppa(e) Mean flux value (mmol/gDCW hr) in vitro in vivo Growth condition Log phase Osmotic stress Stationary phase Sebaceous site T oe web in vitro in vivo Growth condition Log phase Osmotic stress Stationary phase Sebaceous site T oe web Reactions: Reactions: Reactions: Reactions: −5.0 −2.5 0.0 2.5 5.0 7.5 Mean flux value (mmol/gDCW hr) FBA FBA3 L_LACD2 PFK PPCKr PFK_3 TALA 0 2 4 6 8 10 0.25 -0.25 0.00 0.25 0.50 -0.25 0.00 MME MMM2r PPAKr PTA2 PYRCT EX_ppa(e) Mean flux value (mmol/gDCW hr) Mean flux value (mmol/gDCW hr) Cheek Scalp Propionate metabolism Cutibacterium acnes Histidine metabolism Cutibacterium acnes EX_his_L(e) HISt2r HISD IZPN NFLGLNH NFORGLUAH URCN Proline and glutamate metabolism Cheek Scalp EX_pro_L(e) PROt2r PRO1x P5CD .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint 2 7 46 48 2 2 4 6 28 3 4 4 39 7 8 7 2 4 23 13 1 1 2 1 4 1 6 1 3 3 1 1 1 cyclic lactone autoinducer peptide L−lactate oxidase pyruvate oxidase epsilon family phenol−soluble modulin Delta lysin family salivaricin M family lantibiotic plantaricin C family lantibiotic gallidermin/nisin family lantibiotic lacticin 481 family lantibiotic lactococcin 972 family bacteriocin thiazolylpeptide−type bacteriocin thiazolylpeptide−type bacteriocin precursor bacteriocin halocin C8−like domain halocin C8 precursor−like protein Corynebacterium tuberculostearicum Cutibacterium acnes Cutibacterium granulosum Staphylococcus sp. Staphylococcus capitisStaphylococcus caprae Staphylococcus epidermidisStaphylococcus haemolyticus Staphylococcus hominis Staphylococcus pettenkoferi Staphylococcus saccharolyticus Staphylococcus warneri Streptococcus sp.Streptococcus mitisStreptococcus oralis Streptococcus pneumoniae Streptococcus salivarius 10 20 30 40 Number of metatranscriptomes MHS388 MHS349 MHS381 MHS393 MHS598 MHS378 MHS398 MHS376 MHS357 MHS344 MHS421 MHS333 MHS426 MHS407 MHS432 MHS373 MHS442 MHS385 MHS406 MHS414 MHS429 MHS437 MHS440 MHS399 Scalp Cheek Antecubital fossa Volar forearm T oe Web MHS389 MHS390 MHS394 MHS396 MHS379 MHS411 MHS418 MHS422 MHS374 MHS427 MHS384 MHS433 MHS331 MHS443 MHS386 MHS435 MHS328 MHS405 MHS415 MHS430 MHS438 MHS594 MHS439 MHS371 MHS391 MHS341 MHS350 MHS380 MHS355 MHS336 MHS419 MHS423 MHS334 MHS360 MHS408 MHS387 MHS329 MHS416 MHS362 MHS366 MHS441 MHS593 MHS412 MHS382 MHS351 MHS338 MHS337 MHS345 MHS424 MHS335 MHS361 MHS342 MHS364 MHS369 MHS343 MHS365 MHS330 MHS363 MHS368 MHS367 MHS413 MHS392 MHS383 MHS339 MHS377 MHS358 MHS425 MHS375 MHS428 MHS409 MHS434 MHS332 MHS370 MHS410 MHS436 MHS372 MHS347 MHS402 cyclic lactone autoinducer peptide * : p < 0.05 *** : p < 0.001 *** *** * ***Autoinducing peptide Free radical generation Phenol soluble modulinswith anti-microbial potential Bacteriocins andanti-microbial peptides Site-specific enrichment L−lactate oxidase pyruvate oxidase epsilon family phenol−soluble modulin Delta lysin family phenol-soluble modulin γ PSM−delta family phenol−soluble modulin lichenicidin A2 family type 2 lantibiotic salivaricin M family lantibiotic class IIb bacteriocin, lactobin A/cerein 7B family plantaricin C family lantibiotic gallidermin/nisin family lantibiotic lacticin 481 family lantibiotic lactococcin 972 family bacteriocin thiazolylpeptide−type bacteriocin thiazolylpeptide−type bacteriocin precursor bacteriocin halocin C8−like domain halocin C8 precursor−like protein Metatranscriptome detection Figure 5 A TM-score: 0.77894 RMSD: 2.36 Streptomyces papain inhibitor (5ntb_B) Barwin domain containing protein (DNF11_2196) Z-score (DALI): 14.5 Seq. identity: 28.9 % n.s. p = 5.32e−04 n.s. p = 1.48e−02 n.s. 1 2 3 0.2% TritonX 5% BHI S. cap S. epi S. hom S. war 0.2% TritonX 5% BHI S. cap S. epi S. hom S. war pro−IL−1B Relative Light Units n.s. p = 5.76e−04 p = 5.01e−04 n.s. n.s. 0 100 200 300 400 500cleaved IL−1B concentration B C D .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted December 4, 2024. ; https://doi.org/10.1101/2024.12.02.626500doi: bioRxiv preprint