**Background:** The human skin microbiome plays a protective role against pathogens, and bacteriophages are known to modulate bacterial communities and transfer host-specific genes. However, the human skin virome and its role in health, particularly viral-host relationships influencing skin microbiome structure and function, remain poorly understood.
**Methods:** The authors analyzed viral metagenomic data from 60 subjects (42 previously sequenced plus 18 new participants) across three anatomical skin locations (right hand, left hand, scalp) at five time points over six months. Viral particle enrichment was performed using 0.22 µm filtration, followed by whole genome amplification and Illumina HiSeq 250 bp paired-end sequencing. Bioinformatic analysis included MEGAHIT assembly, VIBRANT for bacteriophage annotation and auxiliary metabolic gene (AMG) identification, Kaiju for taxonomic classification, and BLAST against the MEGARes v2.0.0 database for AMR gene identification. Targeted amplicon sequencing of the erm(C) gene (299 bp region) was performed on all 60 subjects (n=894 samples) using Illumina MiSeq 250 bp paired-end sequencing. DADA2 was used for amplicon sequence variant (ASV) identification, and protein variant effect prediction tools (PROVEAN, SIFT, Missense-3D) assessed functional impacts of mutations. Statistical analyses included one-way repeated measures ANOVA, PERMANOVA, and generalized linear mixed models.
**Key Results:** A total of 3,230 bacteriophage contigs >1 kb were identified (3,162 lytic, 68 lysogenic; 485 circular, 2,745 linear). Only approximately 16% of phage contigs could be taxonomically classified, predominantly into the class Caudoviricetes (Myoviridae and Siphoviridae families). Alpha diversity showed no significant difference by gender (P=0.6481) or anatomical location (P=0.0923) within subjects, but was highly significant across subjects (P=1.68×10⁻⁸) and over time within subjects (P=0.0048). Beta diversity (PERMANOVA) showed 15% of variation explained by subject (P=0.001, R²=0.15102) and 14% by anatomical location (P=0.001, R²=0.13965). A total of 648 AMGs were identified, primarily involved in carbohydrate metabolism, amino acid synthesis, and cofactor/vitamin metabolism. Sixteen contigs contained AMR genes meeting strict cutoff criteria (e-value <1×10⁻³, >40% coverage, >80% nucleotide identity, bit score >70), including erm(C), par(E), par(EF), par(C), fus(B), msr(D), msr(E), mef(A), mph(E), rlm(H), mer(B), mer(G), and qac(F). The erm(C) gene was the most abundant AMR gene and showed no significant change in abundance over time (P=0.734) and no significant difference within subjects over time (P=0.2168). Targeted sequencing identified 79 erm(C) ASVs, with ASV_1 having 100% nucleotide and amino acid identity to the reference erm(C) protein (WP_001003263.1). Only functional variants of erm(C) persisted at high abundance across subjects; low-abundance ASVs contained mutations predicted to produce non-functional proteins.
**Clinical Implications:** This study demonstrates that the human skin virome harbors a diverse array of auxiliary metabolic and antimicrobial resistance genes that may enhance bacterial host fitness and persistence. The presence of clinically relevant AMR genes such as erm(C)—which confers macrolide, lincosamide, and streptogramin B resistance and is a leading AMR mechanism in MRSA—within stable phage populations suggests that bacteriophages serve as reservoirs and vectors for AMR gene dissemination on the skin. The finding that only functional erm(C) variants persist at high abundance indicates evolutionary selection for functional AMR genes, potentially contributing to the stabilization of resistant bacterial populations. These results underscore the need for system-based approaches to understand microbiome function and highlight the skin virome as an underappreciated factor in antimicrobial resistance dynamics.