Sixty-five volunteers (33 females, 32 males) with diverse scalp conditions were recruited from Shantou and Putian, China. Participants filled out a questionnaire regarding hair washing frequency, scalp type, and self-assessed condition. Scalp samples were collected aseptically using sterile cotton swabs in conjunction with sterile combs. Each scalp area was thoroughly combed for 5 min, repeated three times, to maximize sample collection. Swabs were immediately immersed in 1 mL DNA/RNA Shield (Jianshibio, China) and stored at −20°C until further processing. The study protocol received ethical approval from Guangdong Technion-Israel Institute of Technology-E20231121001. All participants granted their informed consent before participation enrolment.

Participants were categorized into three groups based on self-assessment questionnaires and visual inspection by a trained dermatologist: (1) Dandruff Oily (DO): Individuals reporting visible dandruff flakes and exhibiting oily scalp conditions. (2) Healthy Oily (HO): Individuals reporting no visible dandruff flakes and exhibiting oily scalp conditions. (3) Healthy Non-Oily (HN): Individuals reporting no visible dandruff flakes and exhibiting non-oily scalp conditions. The exclusion parameters encompassed: any antecedent of cranial dermatoses (such as psoriasis or eczema), utilization of topical or systemic antifungal or antibiotic pharmacological agents within the preceding month, and any additional significant health afflictions that could influence the cutaneous microbiome. All participants supplied documented informed consent before registration. The visual inspection was conducted to objectively identify dandruff flakes.

Genomic DNA was extracted from samples using the ZymoBIOMICS DNA Microprep Kit (Zymo Research, USA) following the manufacturer's instructions, including a 10-min bead-beating step using an MX-C cell disruptor (Dlab, China) to ensure efficient cell lysis.

Amplicons for the bacterial 16S rRNA gene V4 region were generated using the primers GTGCCAGCMGCCGCGGTAA and GGACTACHVGGGTWTCTAAT. For the fungal ITS1 region, the primers used were CTTGGTCATTTAGAGGAAGTAA and GCTGCGTTCTTCATCGATGC. PCR reactions (15 μL final volume) employed Phusion^® High-Fidelity PCR Master Mix (New England Biolabs). PCR products were purified using magnetic beads (Agencourt AMPure, XP Beckman Coulter) and size-selected via gel electrophoresis. Illumina library preparation was performed according to the manufacturer's protocol, with quality assessment by AATI Fragment Analyzer and quantification using Qubit fluorometry and qPCR. Sequencing was performed on an Illumina NovaSeq 6000 platform using a paired-end 250 bp (PE250) protocol.

Raw reads were processed using FASTP (v0.23.1) for quality filtering and adapter trimming (Chen et al., 2018). Paired-end reads were merged using FLASH (v1.2.11) (Magoc and Salzberg, 2011). Amplicon sequence variants (ASVs) were identified using DADA2 within QIIME2 (v2022.2) (https://qiime2.org) (Callahan et al., 2016). The SILVA database was utilized for bacteria taxonomic assignment, while the UNITE database was used for fungi (Yilmaz et al., 2014; Abarenkov et al., 2024). Chimeric sequences were removed using QIIM2′s inbuilt chimera detection algorithm. In addition, to minimize host DNA contamination of the samples, all the metagenomic raw reads were screened and the reads matched to human genome in NCBI database were discarded and removed.

Seven samples per group (HN, HO, DO) were selected from the full cohort (n = 65) for shotgun metagenomic sequencing, matched by age (± 3 years), sex, geographic origin, and scalp condition. DNA from these samples was pooled per group. These seven samples per group were selected to balance statistical power with sequencing costs, as preliminary power analyses suggested this would capture major trends. DNA was fragmented to ~350 bp using a Covaris ultrasonic disruptor. Illumina library preparation and sequencing were performed using a PE150 protocol on the NovaSeq 6000 platform. Raw reads were processed using FASTP (v0.23.1) for quality filtering and adapter trimming. High-quality reads were assembled using MEGAHIT (Li et al., 2015).

High-throughput sequencing of the 16S rRNA gene V4 region yielded an average of 118,498 raw reads per sample (n = 65), while ITS1 region sequencing generated 54,612–154,544 reads per sample (mean = 107,198 reads). The top ten bacterial phyla and thirty most abundant genera are detailed in Supplementary Table S1. Analysis showed Actinobacteriota as the main bacterial phylum in all samples, followed by Firmicutes. Three dominant genera were Cutibacterium, Staphylococcus, and others Lawsonella. Fungal communities were predominantly Basidiomycota (98% average), with Malassezia and Malasseziales_gen_Incertae_sedis as the most prevalent genera.

To investigate microbiome differences across scalp conditions, samples were categorized into three groups: healthy non-oily (HN), healthy oily (HO), and dandruff oily (DO), each comprised of seven samples. Actinobacteriota was more abundant in HN (79.45%) and HO (74.94%) than DO (40.65%, p < 0.0001). Conversely, Firmicutes were more abundant in DO (58.85%) than HN (19.63%) and HO (24.22%, p < 0.0001) (Figure 1A). Other phyla exhibited low relative abundance across all groups. This significant shift in the relative abundance of dominant bacterial phyla suggests a substantial alteration in the microbial ecosystem between healthy and dandruff-affected scalps.

Genus-level analysis (Figure 1B) revealed further distinctions. Staphylococcus was more abundant in the DO group (58.70%) than in HN (19.35%) and HO (23.85%, p < 0.0001). In contrast, Cutibacterium was significantly more abundant in HN (64.98%) and HO (63.24%) groups compared to DO (27.87%, p < 0.0001). Corynebacterium was more abundant in HO (0.43%) and DO (0.41%) than in HN (0.19%). These genus-level variations suggest a potential role for these bacterial taxa in dandruff development or progression.

Fungal analysis (Figure 1C) showed Basidiomycota as the dominant phylum (>98%) in all groups. Ascomycota was more prevalent in HN (1.14%) than in HO (0.18%) and DO (0.08%). Rozellomycota, Mortierellomycota, and Olpidiomycota were present at low abundances across all groups. While significant differences were observed at the fungal genus level (Figure 1D), these variations were less pronounced compared to the observed bacterial shifts. Malassezia represented the majority of fungal communities across all groups (90.14%, 85.05%, and 78.80% in HN, HO, and DO, respectively). The DO group exhibited a higher relative abundance of Malasseziales_gen_Incertae_sedis (21.10%) compared to HN (8.36%) and HO (14.62%).

The composition of the overall taxonomy and abundance patterns resulting from the 16S rRNA and ITS1 amplicon sequencing (Figure 1) largely agreed with findings from shotgun metagenomic sequencing (Supplementary Figure S1). Both methodologies indicated a substantial increase in Staphylococcus and a decrease in Cutibacterium abundance within the DO group compared to the healthy cohort.

Alpha diversity analysis assessed within-group diversity using Chao1, observed operational taxonomic units (OTUs), dominance, Good's coverage, Pielou's evenness, Simpson's, and Shannon diversity indices (Supplementary Table S1). No significant differences in alpha diversity existed among the three groups (HN, HO, DO), indicating similar species richness and evenness for bacterial and fungal communities (Figures 2A, C).

Beta diversity analysis using Bray-Curtis dissimilarity compared microbial community structures across groups. Bacterial community structure differences were observed between groups (p < 0.05, Figure 2B). This contrasted with the fungal communities, where no statistically significant differences were found across groups (Figure 2D).

Principal coordinate analysis (PCoA) of Bray-Curtis dissimilarities from 16S rRNA data (Figure 2E) showed a clear separation between healthy groups (HN, HO) and the dandruff oily (DO) group. Analysis of similarities (ANOSIM) confirmed significant differences in bacterial community structure between DO and both HN (p = 0.015) and HO (p = 0.005). This contrasted with the fungal communities, where no statistically significant differences were found across groups (Figure 2D).

Pairwise comparisons of bacterial genera was conducted through genus-level t-tests (Figures 3A, B). The results were further illustrated using volcano plots (Figures 3C, D), which combined p-values with fold changes in abundance. A noteworthy difference in the average abundance of Cutibacterium and Staphylococcus was found between the HO and DO groups (p < 0.001 for both genera). Cutibacterium was significantly more abundant in the HO group, while Staphylococcus predominated in the DO group. The volcano plot (Figure 3C) confirmed this, showing a substantial positive fold change [–log10(p-adj) > 4] for Cutibacterium in HO relative to DO, indicating strong upregulation in healthy oily scalps. Conversely, Staphylococcus showed a significant negative fold change in HO compared to DO. Comparative analyses of the DO and HN groups (Figure 3B) indicated a notably higher prevalence of Staphylococcus in the DO group (Figure 3D). While Cutibacterium showed numerically lower abundance in DO compared to HN, the difference was not statistically significant.

Functional profiles were predicted using PICRUSt2, which inferred pathways from 16S rRNA gene sequencing and the KEGG database. Metabolic pathways made up about 50% of the predicted functions: 50.87% (HN), 50.72% (HO), and 49.84% (DO). Genetic Information Processing constituted the second most prevalent category, comprising 18.65%, 18.63%, and 18.15% for the HN, HO, and DO groups, respectively. This was succeeded by Environmental Information Processing, which accounted for 15.91%, 16.01%, and 15.81% for HN, HO, and DO, correspondingly. These three functional categories represent the primary areas of predicted functional enrichment in the microbial communities across the three scalp groups (Figure 4A). In our analysis, the functional predictions from DIAMOND alignment (Supplementary Table S1) and PICRUSt2 were broadly consistent and both methods showed similar trends.

To examine inter-genus relationships, co-occurrence networks at the genus level were established for each respective genus group. Distinct interaction patterns emerged across groups. Within the HN group, the predominant genera, namely Cutibacterium, Staphylococcus, and Lawsonella, demonstrated intricate interrelationships with various other genera, as depicted in Figure 4B. A negative correlation between Cutibacterium and Lawsonella suggests potential competitive or inhibitory interactions. The HO group displayed similar interaction patterns but also included Pseudonocardia within its core network, and a negative correlation was observed between Lawsonella and Pseudomonas, a notable difference from the HN group (Figure 4C). In the DO group, both Cutibacterium and Lawsonella showed individual negative correlations with other genera (Methylobacterium-Methylorubrum and Porphyromonas, respectively). Pseudomonas exhibited extensive connections with various genera, including Streptococcus, highlighting its significant ecological role. Positive correlations were observed between Staphylococcus and Pseudonocardia (p = 0.005), as well as Acinetobacter (p = 0.029). The absence of correlations between Staphylococcus and Cutibacterium or Lawsonella in the DO group suggests niche differentiation or specialization (Figure 4D).

Shotgun metagenomic sequencing revealed notable differences in the scalp microbiome composition between healthy individuals and those affected by dandruff individuals. Taxonomic annotation, visualized using Krona (Supplementary Figure S1), provided a high-level overview of taxonomic distribution at the phylum and class levels. The dandruff oily cohort showed a bacterial composition comprising 65% of the microbiome, with increased Staphylococcus and decreased Cutibacterium abundance (Figure 5A).

Species-level analysis (Figure 5B) highlighted specific taxa differentiating the groups. Cutibacterium acnes showed a markedly higher prevalence in the HN (7.21%) and HO (5.45%) groups compared to the DO group (0.78%). Similarly, Cutibacterium namnetense showed higher relative abundance in HN (1.13%) and HO (0.89%) compared to DO (0.36%). In contrast, the DO group showed a higher prevalence of Staphylococcus aureus at 1.31%, compared to HN and HO. Malassezia restricta, frequently associated with dandruff, exhibited higher abundance in the HO (4.32%) and DO (3.39%) groups relative to HN (2.78%). Lactobacillus crispatus, recognized for its probiotic properties, was more prevalent in the DO group (0.83%) compared to HN and HO. These variations in species abundance may significantly influence alterations in skin barrier function observed in individuals with dandruff.

Further analysis involved a statistical heatmap of unigene annotation numbers and a species abundance clustering heatmap (Figures 5C, D) to visualize differences in microbial composition across the samples. Malassezia restricta exhibited consistently high unigene abundance across all groups (2,497, 2,526, and 2,510 in HN, HO, and DO, respectively). Conversely, Cutibacterium acnes abundance was significantly higher in HN (887) and HO (893) groups compared to DO (565). The species abundance heatmap (Figure 5D) indicated that the DO group exhibited greater abundances of Staphylococcus, Pseudomonas, and Cutibacterium granulosum. In contrast, Cutibacterium acnes and Cutibacterium namnetense were more prevalent in the HN group compared to the HO and DO groups. Malassezia restricta also showed greater relative abundance in HO compared to HN and DO.

Shotgun metagenomic sequencing data allowed for the investigation of the functional potential of the scalp microbiome across the three groups (HN, HO, DO). Level 1 KEGG pathway analysis, consistent with previous PICRUSt2 predictions (Supplementary Figure S2), revealed distinct functional profiles. Bray-Curtis hierarchical clustering (Figure 6A) indicated that HN and HO samples were functionally similar. Compared to DO, HN, and HO groups showed higher relative abundance of carbohydrate/amino acid metabolism, nucleotide/energy metabolism, and membrane transport pathways (though statistical testing was precluded by sample pooling). Furthermore, HN and HO groups showed greater representation of nucleotide and energy metabolism pathways compared to the DO group. The “Environmental Information Processing: Membrane Transport” KEGG category was significantly enriched in both HN and HO groups, suggesting differences in nutrient acquisition and environmental responsiveness between these groups and the DO group. Similarly, the “Genetic Information Processing: Translation” pathway was significantly enriched in HN and HO groups.

To further refine the functional analysis, Enzyme Commission (EC) numbers were annotated to identify specific enzymes. Hierarchical clustering of EC numbers (Figure 6B) showed that EC number 7.6.2.4, linked to ABC-type fatty-acyl-CoA transporters, was mainly found in the DO group, indicating a possible role in fatty acid transport. In contrast, EC numbers 3.6.4.13 (an ATPase linked to energy metabolism and transport) and 3.1.3.16 (a phospholipase connected to lipid metabolism and membrane dynamics) were significantly more common in the HO group.

Orthologous gene clusters were analyzed across scalp groups using the EggNOG database to assess functional differences. Hierarchical clustering (Figure 7A) showed that the DO group displayed increased activity in functions associated with replication, recombination, and DNA repair, indicating improved genomic maintenance and stability. In contrast, the HN group showed greater activity in transport and metabolism pathways, particularly membrane biogenesis, reflecting adaptations for nutrient uptake and cellular integrity. The HO group showed notable activity in post-translational modification, protein turnover, chaperone functions, secretion, vesicular transport, and cell cycle regulation. These variations highlight the distinct ecological roles and adaptive strategies of microbial communities within each scalp group.

Further analysis using level 2 hierarchical clustering (Figure 7B) revealed more specific functional enrichments. The DO group showed significant enrichment in transposition-related functions, with abundant transposases and transposable element-associated domains (L1, LINE1), indicating a high degree of genomic plasticity. The HN group displayed significant enrichment in carbohydrate transport functions, primarily via ABC transporters and the major facilitator superfamily, suggesting efficient carbohydrate metabolism and utilization. The HO group exhibited a more balanced functional profile, with moderate representation of both transposition-related and carbohydrate transport functions, along with rRNA binding and transcriptional regulation, consistent with a more generalist metabolic strategy compared to the HN and DO groups.

Hierarchical clustering of CAZyme family abundances (Figure 8A) revealed distinct functional adaptations across the three groups, potentially correlating with sebum composition and dandruff pathogenesis. The DO group showed a relatively lower abundance of polysaccharide lyases (PLs), suggesting reduced capacity for complex polysaccharide degradation. The HN group showed more glycosyl hydrolases (GHs), especially for degrading simple sugars from host-derived (skin cells) and environmental sources (e.g., skincare products). Increased levels of carbohydrate esterases (CEs) in the HN group, and comparable levels of glycosyl transferases (GTs) in both HN and HO groups, suggest enhanced capabilities for utilizing lipid-associated carbohydrates and modifying host-produced lipids. The HO group showed increased relative abundance of auxiliary activities (AAs) and PLs, indicating a greater capacity for complex polysaccharide degradation, which may contribute to the breakdown of stratum corneum components, potentially mitigating barrier disruption and inflammation associated with dandruff.

Hierarchical clustering utilizing Bray-Curtis dissimilarity (Figure 8B) showed that the healthy groups (HN and HO) exhibited higher community cohesion and resilience compared to the DO group. Although essential virulence factors (VF0317, VF0842, and VF0302) existed in every group, they were found to be more plentiful in the healthy groups. The higher abundance of VF0842 in the DO group suggests a specific role in dandruff pathogenesis. Other virulence factors (VF0738, VF0739, VF0758, VF0755, VF0298, VF0849, and VF0858) were present across groups but exhibited lower prevalence in the DO group.

Analysis of pathogen-host interactions using the PHI database (Figures 8C, D) revealed a higher relative abundance of factors associated with skin infections in the DO group, indicating a skewed microbial community potentially contributing to dandruff. At level 2, specific pathogen-host interactions were highlighted. The much higher levels of Staphylococcus aureus in the DO group indicate a connection to inflammation and infection, upset the balance of the skin microbiota. Conversely, the HO group exhibited increased levels of fungi (Ustilago maydis and Candida albicans), which may contribute to scalp health through competition with pathogenic organisms. The HN group exhibited a more balanced microbial profile, indicative of a healthier scalp environment.

Figures 9A–D show significant differences in the relative abundance of ARGs and mobile genetic elements (MGEs) among the three scalp samples groups. Figure 9A shows the relative abundance of key ARGs. The vanT gene (vanG cluster) consistently exhibits high relative abundance, particularly in healthy groups (76.41% in HO and 86.16% in HN of the total ARG pool), suggesting widespread potential for vancomycin resistance. Other ARGs (sdrM, norC, mdeA, fosBx1, sepA, mgrA, and inuA) account for the remaining 15–25% of the ARG pool. While these genes maintain similar proportions in the DO group (8.23–15.6%), with vanT, sdrM, and norC each representing ~15.6%, many are associated with Staphylococcus aureus, suggesting potential contributions to resistance mechanisms.

Figure 9B depicts the relative abundance of integrons. The integron CP012744.1 constitutes ~40–50% of the total integron pool, indicating its significant role in the microbiome's genetic landscape. This integron is significantly more abundant in the DO group (53.81%) compared to HN (40.01%) and HO (41.45%), suggesting an adaptive response to selective pressures associated with dandruff. In contrast, integron LT996886 is nearly five times more prevalent in the healthy groups, possibly reflecting a protective role in maintaining scalp health. The diverse integron composition further underscores the complexity of microbial interactions and resistance mechanisms.

Figure 9C shows the relative abundance of the top 10 insertion sequence (IS) elements. In HN samples, ISArsp12 (68.01%) and ISSau6 are most prevalent. ISSau6, linked to Staphylococcus aureus, is significantly more abundant in DO (54.71%) compared to HN (21.28%) and HO (17.67%), suggesting a microbial shift during dandruff development that may contribute to inflammation. ISCge1, primarily found in HO (45.65%), exhibits lower frequency in HN (0.36%) and DO (15.12%), suggesting a potential protective role. Other IS elements show subtle variations.

Figure 9D shows the relative abundance of plasmids. Plasmid NZ_CP053353.1, originating from Staphylococcus aureus, is significantly more abundant in DO (39.49%) compared to HN (9.24%) and HO (11.77%), suggesting a link to dandruff-associated microbial dynamics. Healthy groups show a relatively even distribution of plasmids NC_021086.1, NZ_CP017041.1, and NZ_LN868939.1 (36–40% of total plasmid pool), potentially indicating a more stable plasmid composition reflecting a balanced ecosystem. The remaining plasmids exhibit considerable diversity and variation across groups.