Evolutionary and functional implications of hypervariable loci within the skin virome

Analysis details and availability

All associated source code and explanatory README files are available for review at the following GitHub repository: https://github.com/Microbiology/ViromeVarScripts.

Data acquisition and quality control

The primary skin virome dataset was acquired from SRA accession: SRP049645 (Hannigan et al., 2015) and includes sequences from samples collected at the second and third time points. The secondary dataset was obtained from Oh et al. (2014) (SRA BioProject: 46333). Retroauricular crease samples were downloaded from the NCBI SRA BioProject: 46333. For samples from both the primary and secondary dataset, sequences were trimmed with the FASTX-Toolkit (version 0.0.14), using a quality score cutoff of 33. Remaining sequences with similarity to the human genome were removed using the standalone DeconSeq toolkit (version 0.4.3) (Schmieder & Edwards, 2011).

Contig assembly and taxonomic identification

Contigs from the primary dataset were obtained from the published Figshare source (DOI 10.6084/m9.figshare.1281248). Contigs from both the primary and secondary datasets were separately assembled using the Ray metagenomic assembly software, specifying a minimum contig length of 500 bp and otherwise default parameters (v2.3.1) (Boisvert et al., 2012). Within each dataset, sequences from all samples were combined prior to assembly to facilitate the most complete contig assembly. Contig coverage was determined by aligning sequences back to the contigs with the bowtie2 toolkit (v2.1.0; seed substring length of 25 and one mismatch allowed in alignment) (Langmead & Salzberg, 2012). Quantification of reads mapping back to contigs was obtained by parsing bowtie2 output using Perl and BASH scripts presented in the supplemental source code. Coverage was calculated using the number of reads mapping to each contig. The blastn program from the NCBI Blast+ toolkit (version 2.2.0) was used to determine similarity of contigs to virus reference genomes (Camacho et al., 2009). Contigs were blasted against a previously described virus-specific genome reference database, which is a subset of the EMBL reference genome database (UniProt Consortium, 2014; Hannigan et al., 2015). A similarity threshold of e-value < 10⁻³ was used, and sequences with multiple potential identities were resolved by using only hits with the lowest e-values. Although this was the minimum threshold, the contigs of interest exhibited e-values < 10⁻³.

Phylogenetic analysis

We constructed phylogenies using the L1 capsid gene for HPV (Ma et al., 2014) and the large terminase subunit for the Staphylococcus and Propionibacterium phages (Gutiérrez et al., 2013; Ma et al., 2014) as phylogenetic marker genes. For reference, we used the PAVE reference L1 genes (https://pave.niaid.nih.gov/, accessed 2015-06-03) (Van Doorslaer et al., 2013). The large terminase subunit references for Staphylococcus and Propionibacterium phages were from the NCBI gene sequence database (Staphylococcus phage: accessed 2015-09-14, search term: ((phage terminase large subunit staphylococcus)) AND “viruses”[porgn:__txid10239] NOT “ORF” NOT “hypothetical” NOT “putative;” Propionibacterium phage: accessed 2015-09-15, search term: ((phage terminase large subunit propionibacterium)) AND “viruses”[porgn:__txid10239]). To extract the phylogenetic marker genes from the virome contigs, we determined which open reading frames (ORFs) matched the reference genes by nucleotide similarity (nucleotide blast, e-value 1e−10). Only ORFs longer than 1.2 kb were included in the analysis. The average reference gene lengths were all longer than this threshold (average reference gene length: HPV = 2,519 bp, Staphylococcus phage = 1,307 bp, Propionibacterium phage = 1,511 bp). Contig and reference marker genes were aligned using the Smith–Waterman algorithm and 1,000 iterations as implemented by the mafft aligner (v7.221) (Katoh & Standley, 2013). Phylogeny was constructed using RAxML (version 8.1.21) (Stamatakis, 2014). The phylogenetic tree was visualized using Figtree (Rambaut, 2006).

Identification of temperate phage contigs

As has been described previously, we identified temperate (lysogenic) phage contigs using three genomic markers: contig nucleotide similarity to (1) phage integrase genes, (2) prophage genes within the ACLAME prophage database, and (3) bacterial reference genomes. We performed a blastx alignment (e-value 1e−10, percent identity threshold 75%) of the genes within the ACLAME prophage database (Leplae, Lima-Mendez & Toussaint, 2010), a blastx alignment with integrase genes from Uniprot database, and a blastn alignment of the Staphylococcus phage contigs to Staphylococcus bacterial reference genomes. Integrase genes were obtained from the online Uniprot database using the search term “organism:phage gene:int NOT putative.” Staphylococcus reference genomes were obtained from the NCBI nucleotide database using the search term “‘Staphylococcus’[Organism] AND ‘complete genome’[Name] NOT phage[All Fields] NOT contig[All Fields] NOT (‘unidentified plasmid’[Organism] OR plasmid[All Fields]) AND (bacteria[filter] AND biomol_genomic[PROP]).” Both were accessed December 22, 2016. Together this allowed us to detect regions of contigs that demonstrated a high similarity to temperate phage gene signatures.

Identification of hypervariable loci

The bowtie2 alignments of reads to viral contigs were formatted (e.g., conversion from binary to ASCII format) and then single nucleotide polymorphisms (SNPs) were called using VarScan (v2.3.7) (Li et al., 2009; Koboldt et al., 2012). The “pileup2snp” program from VarScan was used with a minimum minor allele frequency threshold of 1%, a read depth of 8, and a minimum of two supporting reads for variant calls. Indels were excluded.

To identify hypervariable loci, we used a geometric distribution based statistic approach as described previously (Zheng et al., 2010), which, compared to sliding window searches and other similar methods, has the advantage of avoiding boundary difficulties and variations within contigs. We used a geometric distribution to model the probability of achieving two SNPs separated by a specified non-SNP nucleotide distance. Each between-SNP distance was associated with a probability and the probability of a particular distance occurring by randomly sampling was less than 5%. Thus, we identified a range of SNP distances as significantly less than background if they occurred within our dataset less than 5% of the time.

Protein family domain identification within hypervariable loci ORFs

Protein family domains were identified in ORFs that contained hypervariable loci. The subset of translated virus ORFs that contained hypervariable loci were aligned to the standard Pfam protein family domain database using hmmscan within the HMMer toolkit (version 3.1) and GA gathering bit score thresholds (Finn, Clements & Eddy, 2011).

Prediction of single amino acid variant effect on phenotype

The SuSPect algorithm was employed to predict the likelihood of SNP-associated single amino acid variants (SAVs) impacting phenotype. We used SuSPect to create a matrix of likelihood scores for every possible SAV at every position in the ORFs that contained hypervariable loci. This matrix was used as a reference to quantify the likelihood of each hypervariable loci SNP to impact the resulting phenotype. The significance of the score differences between viruses was calculated using a Wilcoxon rank-sum test.

Evolutionary pressure of hypervariable loci and virus genomes

We assessed the evolutionary pressure of a gene using the pN/pS ratio as in Formula 1, where M_N and M_S represent the observed number of non-synonymous and synonymous SNPs, respectively. These values were normalized by the total number of possible non-synonymous or synonymous substitutions (N_i and S_i, respectively), in order to avoid potential codon usage bias. Furthermore, to normalize for sequence coverage of the SNPs and prevent extreme values, a pseudocount value of an arbitrarily small number was added to the M_N and M_S values, which was defined as half of the square root of the median sequence coverage of SNPs within the region of interest (C_M). The pseudocount approach was used to prevent infinite and illegal values when M_N or M_S had zero values, thus allowing consideration of otherwise infinite or ignored data points. For example, an absence of synonymous mutations would result in an infinitely large value (dividing by zero) thus forcing exclusion of the data point. Our approach preserves this data and allows us to draw conclusions from the largest possible dataset, and has been shown to be effective in previous studies (Novaes et al., 2008; Bajgain et al., 2011).

The formula used to calculate the pN/pS ratio for a gene. M_N is the number of non-synonymous SNPs within the gene and M_S is the number of synonymous SNPs found within the gene. Each mutation value is normalized for the likelihood that the result would have happened by chance, calculated as the sum of the proportions of nucleotides that would have resulted in either a non-synonymous or synonymous mutation. To calculate this proportion, the possible non-synonymous mutations (N) and synonymous mutations (S) at position i within the gene are expressed as a fraction of the three possible alternate nucleotides. SNP counts were smoothed as pseudo-counts using the median SNP sequence coverage (C_M).

This analysis is similar to the dN/dS calculations often performed to estimate degrees of natural selection among genomes (Nishida, Frith & Nakai, 2009; Schloissnig et al., 2013). It is important to note however that such an analytical approach would be inappropriate for this type of sample set because the nucleotide variations are not assignable to isolated strains, which prevents haplotype identification that is a necessary component of dN/dS calculations. The pN/pS calculation does not assume haplotypes, and is therefore appropriate for metagenomic datasets.

To estimate the selective pressure on hypervariable loci, the locations of the hypervariable loci were extracted, along with their immediately adjacent regions, using a Perl script as presented in the supplemental source code. Adjacent regions are defined as the genomic regions that are two times the length of the hypervariable loci and located immediately before and after the hypervariable loci. These positions were used with the contig sequences, SNP call data, and pN/pS calculator to estimate their selective pressure. Hypervariable regions outside of coding regions were not considered.

Calculation of the overall selective pressure on virus contigs was performed in a similar approach to the hypervariable loci selective pressure. Predicted ORFs were first extracted from the contigs using the Glimmer3 toolkit (v3.02) (Delcher et al., 2007). The predicted ORFs, along with the contig sequences, SNP profile, and pN/pS calculator were used to calculate the overall selective pressure on each gene within each contig. The distributions of selective pressures observed for each gene were observed as categorized by virus type.

Amino acid frequency, charge, and polarity

Amino acid abundance profiles were calculated while correcting for the random probability of that substitution. More specifically, each value was weighted for the number of nucleotides that result in the same amino acid as weighted value = ((number of nucleotide substitutions resulting in same amino acid)/3)⁻¹. Relative abundance was calculated as the sum of the corrected frequencies. Charge and polarity were determined using a simple table of known amino acid properties. Differences in profiles between viruses were calculated using a chi-square test.

Diversity generating retroelement identification

We identified potential DGRs by collecting assembled contigs that contained ORFs similar to known reverse transcriptase genes, and a duplicated nucleotide region less than 150 bp in length. Reverse transcriptase ORFs were identified using blastx (e-value < 10⁻⁵) and the Uniprot reference reverse transcriptase sequences (http://www.uniprot.org/uniprot/?sort=score&desc=&compress=yes&query=%22reverse%20transcriptase%22%20(phage%20OR%20virus)&fil=&format=fasta&force=yes). Repeat regions were identified by comparing each contig to itself with tblastx (e-value < 10⁻⁵⁰) and were filtered using custom scripts to remove duplicates and regions longer than 150 bp. DGR candidates were removed if they contained no hypervariable loci or if the variable region was not within a predicted ORF.

Diversity generating retroelements were visualized in the Integrated Genomic Viewer using the DGR cassettes and bowtie2 aligned sequences described above. The linkage disequilibrium was calculated using a custom Perl script for formatting and the “LDheatmap” and “genetics” R packages for analysis and visualization (Shin et al., 2006; Warnes et al., 2003). The linkage disequilibrium for each pair of SNPs was calculated as the squared allelic correlation (R²).

Comparison of primary analysis to validation dataset

Near identical contigs were identified between the primary and secondary validation dataset by aligning the two individually assembled contigs to each other using bowtie2, with a specified seed length of 25 and up to one seed mismatch. Sequences from our primary dataset and the Oh et al. (2014) dataset were aligned to the near identical contigs. These alignments were used to identify shared SNP locations between our dataset and the Oh et al. (2014) dataset. We quantified shared SNP location as percent of our primary analysis SNPs whose location was identical to those of SNPs in the secondary validation dataset. As a control, we compared these results to a simulated dataset where SNP position was randomly assigned. The example SNP alignment over the circular contig was generated using Genious (Kearse et al., 2012).

Diversity of skin viruses

We evaluated genomic variability associated with dsDNA skin viruses using a previously published human skin virome metagenomic dataset, consisting of 260,714,906 high quality sequences assembled into >76,000 contigs from 16 individuals (SRA accession: SRP049645) (Hannigan et al., 2015). We relied on database virus annotation to identify the taxonomic groups whose contigs had the overall highest confidence matches to reference genomes. Because greater sequencing coverage allows for more refined detection of variable nucleotides (Schloissnig et al., 2013), we focused our analysis on taxa whose de novo assembled contigs had sufficient coverage (greater than 10×). Contigs meeting these criteria were identified as Propionibacterium phages (contig count = 45), Staphylococcus phages (contig count = 319), and human papillomaviruses (HPVs; contig count = 56; Fig. 1A), representing a total of 420 contigs to be used out of the more than 76,000 total contigs in the study (Fig. S1). All of the contigs from these three taxa were used in our analysis, including those below our coverage threshold, since contigs can have regions of high coverage despite an average low coverage. More specific, secondary filtering was done while identifying SNPs. Some contigs were identified as Pseudomonas phages or Enterobacteria phages, but these taxa excluded from the analysis because their annotations were lower confidence and contig representation was minimal (Fig. 1A). The uneven virus population coverage is potentially a reflection of our inability to taxonomically identify the majority of virome sequences, and as a result we lose this information to the viral “dark matter” (Hannigan et al., 2015).

We evaluated the diversity of the skin viruses by constructing phylogenic trees based on conserved viral genes. Only fully assembled ORFs were considered of >1.5 kb for HPV, and >0.15 kb for Staphylococcus and Propionibacterium phages. Similar to previous studies, we used the L1 major capsid gene to classify HPV strains (de Villiers et al., 2004; Ma et al., 2014). The terminase large subunit gene was extracted from Staphylococcus phage contigs to construct phylogeny as described previously (Gutiérrez et al., 2013; Ma et al., 2014). Because this gene is used for phylogeny of a variety of phages, we attempted to construct Propionibacterium phage phylogeny in a similar manner (Ganz et al., 2014; Li et al., 2014), but were ultimately unsuccessful due to the lack of a full-length de novo assembled reference genes in the dataset.

Most skin HPVs were identified as gamma HPVs, the prototypical cutaneous HPV class (Fig. 1B) (Mistry, Wibom & Evander, 2008). Few contigs were identified as beta and Mu/Nu HPVs, and none were identified as alpha HPVs. This is consistent with data from the Human Microbiome Project cohort (Ma et al., 2014).

Fewer Staphylococcus phage marker genes were identified, compared to HPVs, likely because Staphylococcus phage genomes are orders of magnitude longer than HPV genomes, thereby decreasing the probability that contigs covered the entire genome. Because multiple displacement amplification (MDA) was not used to create this dataset, there is no MDA-associated bias toward small circular genomes. The Staphylococcus phage contigs belonged to the Phietalikevirus genus and orphan virus groups (those that have not yet been classified) (Fig. 1C). Of the Staphylococcus phage contigs identified, 49.6% (123 out of 248 contigs) were predicted to be lysogenic, based on similarity to lysogenic phages in the ACLAME database, integrase genes in the Uniprot database, and Staphylococcus reference genomes from the NCBI nucleotide database, as described previously (Leplae, Lima-Mendez & Toussaint, 2010; Minot et al., 2011; Hannigan et al., 2015). This is a minimum estimate of contig lysogeny, as some of the other contigs may have lysogenic signatures that we failed to identify. Furthermore, because this classification strategy is based on blast assignments, it may result in false positives if genes in the database are homologous to genes present in lytic phages.

Hypervariable loci within the skin virome

We implemented a geometric distribution-based approach to identify regions of high genomic diversity, as in (Zheng et al., 2010). Regions within each contig that contained a significantly higher frequency of SNPs over the stochastic background were identified as viral hypervariable loci. Significance was defined as the frequency of SNPs having less than a 5% chance of randomly occurring, given the geometric distribution of the dataset. HPVs and Staphylococcus phages maintained 106 and 465 hypervariable loci, respectively. We were unable to detect hypervariable loci in the Propionibacterium phage population.

To determine the virus protein family domains hosting hypervariable loci, we used the hidden Markov model analysis implemented by HMMer (Finn, Clements & Eddy, 2011). Hypervariable loci-containing HPV genes include E6, E2, and E1 genes, which are associated with infectious gene expression, and the L1 major capsid gene, which is involved in tropism and host immune evasion (Table S1). The L1 major capsid protein is also a target in contemporary, widely used HPV vaccines (Schiller & Lowy, 2012). Hypervariable loci were detected in a variety of Staphylococcus phage genes with predicted functions related to tropism, host immune evasion, and utilization of host resources (Table S2).

Selective pressures on hypervariable loci

We evaluated the selective pressures on virus genes by calculating the pN/pS ratio of non-synonymous SNPs (pN) to synonymous SNPs (pS) within each virus taxa (Schloissnig et al., 2013). This was used as an alternative to dN/dS because dN/dS assumes haplotype information which cannot be fulfilled by metagenomic data (Schloissnig et al., 2013). In the pN/pS calculation, neutral evolution is defined as an equal frequency of synonymous and non-synonymous polymorphisms. Selective pressure favors non-synonymous mutations, resulting in increased pN/pS ratios. Purifying selection has the opposite effect. Because the existing model (Schloissnig et al., 2013) is susceptible to stochastic effects and extreme outliers (e.g., division by zero when pS = 0), we added a pseudocount correction (Formula 1).

We determined whether hypervariable loci are in fact loci of focused selective pressure by comparing pN/pS values of the loci to the adjacent genomic regions. pN/pS values of hypervariable loci were significantly lower than adjacent regions in both HPV (median: adjacent = 1.0, hypervariable loci = 0.21; p-value = 3.4e−17) and Staphylococcus phage (median: adjacent = 0.76, hypervariable loci = 0.41; p-value = 1.8e−40) genomes, suggesting purifying selection and a propensity to maintain existing protein sequences (Figs. 1D). HPV hypervariable loci were under significantly greater purifying selection than those of Staphylococcus phages (median: HPV = 0.21, Staphylococcus phage = 0.41; p-value = 4.64e−9) (Fig. S1). Furthermore, not only are the pN/pS values of the hypervariable significantly lower than their adjacent regions, but very few of the loci have a pN/pS value greater than one.

To evaluate whether the observed selective pressure in HPV and Staphylococcus virus communities is genome-wide or localized to hypervariable loci, we quantified the selective pressure on each virus’ genome by calculating the overall pN/pS ratio including hypervariable loci and non-hypervariable loci SNPs. We observed nearly neutral pressure across HPVs and Staphylococcus phages that mirrored pressures to those observed in the regions adjacent to the hypervariable loci (median: HPV = 1.0, Staphylococcus phage = 0.81, p-value = 3.2e−5) (Fig. S2).

Functional implications of targeted substitutions within hypervariable loci

In order to evaluate the specific nucleotide changes occurring at hypervariable loci, as well as to evaluate the implications of specific nucleotide polymorphisms, we quantified the frequency of individual nucleotide substitutions within hypervariable loci. A>C and T>C substitutions were most frequent in HPV hypervariable loci (Fig. 2A). Staphylococcus phages exhibited a significantly different substitution profile (p-value = 0.00018, chi-square test), with the most common substitutions being A>G and G>A transitions (Fig. 2B). HPV and Staphylococcus phage substitutions were more likely to be transitions, with a transition/transversion (ti/tv) ratio of 3.25 and 2.02, respectively.

We predicted how hypervariable loci SNPs might affect protein functionality by evaluating patterns of the amino acid substitutions while correcting for the random chance that the substitution will occur. The most frequent non-synonymous amino acid substitution in HPVs was glycine (consensus amino acid) to valine (variant amino acid, Fig. 2C). While these amino acids are (non-polar and hydrophobic), glycine is less hydrophobic than valine. The most frequent non-synonymous amino acid substitution in Staphylococcus phages was proline to leucine (Fig. 2D), a substitution between a non-polar cyclic amino acid and an aliphatic straight chain amino acid. Profiles of amino acid substitution were significantly different between HPVs and Staphylococcus phages (p-value = 0.0021; chi-square test).

Amino acid polarity and charge were largely maintained in HPV hypervariable loci (Figs. 2E and 2G). In instances of altered charge, visual inspection suggests the most frequent changes were from neutral to positive or negative charge, or positive to neutral charge. Consensus acidic polar residues were not associated with polymorphisms. Staphylococcus phage community hypervariable loci appeared to be under weaker substitution selection, with a greater diversity in amino acid charge and polarity (Figs. 2F and 2H) compared to HPV. Patterns of substitution charge and polarity were not significant (p-value > 0.5; chi-square test) when comparing the entire HPV to Staphylococcus phage substitution profiles.

We reinforced the observed functional implications of hypervariable loci by predicting the effects of their associated SAVs on gene phenotype using the support vector machine algorithm implemented in SuSPect (Yates et al., 2014). This method assigns a deleterious score to each hypervariable loci SNP-associated SAV, with 0 representing a neutral SAV and 100 representing a SAV with high likelihood to impact phenotype. These scores are based on the predicted impact of the SAV on the tertiary and secondary structure of the resulting protein, the location of the SAV within the resulting protein (surface vs core), and whether the SAV has previously been associated with altered protein–protein interactions. Both Staphylococcus phages and HPVs have an abundance of SNPs associated with SAVs predicted to be deleterious (deleterious scores approaching 100) (Fig. 3). The HPV SNPs were predicted to be significantly more likely to impact phenotype than the Staphylococcus phage SNPs (median: HPV = 45, Staphylococcus phage = 17; p-value < 2.2e−16), suggesting that SNPs impact functionality differently between viruses.

Diversity generating retroelements as a mechanism for targeted hypervariability

Diversity generating retroelements are a genetic system used by bacteriophages (as well as bacteria and archea) to promote targeted hypervariability in genes (Doulatov et al., 2004). While DGRs are complex and consist of many components, at their most basic they can be identified as elements consisting of a reverse transcriptase gene and a repeated nucleotide sequence of length <150 bp that is found in two separate locations of the genome (Doulatov et al., 2004; Minot et al., 2012; Schillinger et al., 2012), termed the template region and the variable region. The template region is transcribed, then reverse transcribed in an error-prone fashion. The resulting cDNA is then integrated into the variable region, introducing base substitutions. Targeted hypervariation impacts functions including broadened host cell tropism by mutagenizing a phage tail fiber gene (Doulatov et al., 2004).

We thus sought to identify candidate DGR cassettes within our viral contigs. We defined the candidate cassettes as pairs of non-overlapping regions with similar nucleotide sequences (tblastx of contigs against themselves, e-value < 10e−50) and co-localized on a contig containing a predicted virus/phage reverse transcriptase gene. We only considered cassettes that were located within a predicted viral gene, contained at least one hypervariable locus in their variable region, and exhibited truly random variation (different between reads). Based on these criteria, we identified one Staphylococcus phage DGR candidate that contained hypervariable loci. We also identified five other DGR candidates that were associated with hypervariable loci outside of predicted genes or that failed to demonstrate linkage disequilibrium, suggesting an association with cryptic genes or pseudogenes.

For the Staphylococcus phage DGR candidate with hypervariable loci, we calculated the linkage disequilibrium associated with the variable nucleotide positions to infer whether the DGR was active or inactive (e.g., an evolutionary artifact). The DGR cassette had unlinked nucleotide variation, which was supported by low levels of linkage disequilibrium (squared allelic correlation R²) between SNP pairs in the variable region (Fig. 4). In this cassette, the template region has less frequent blocks of linkage equilibrium (unlinked variants) while the variable region was associated with greater linkage equilibrium. Together this suggests the observed Staphylococcus phage is active. The variable region was associated with a gene of unknown function.

Skin virome variability patterns and SNP locations are reproducible across different datasets

We repeated our analyses in a separate, independently collected dataset from another research group (SRA BioProject: 46333) (Oh et al., 2014) to determine the generalizability of our findings. We analyzed metagenomic sequence data of skin specimens that were collected from the retroauricular crease without initial purification of virus like particles. Consistent with our primary analysis, Staphylococcus and Propionibacterium phages were identified as having the highest coverage and similarity to reference genomes (Fig. 5A). Pseudomonas phages were identified but were in the minority and had low coverage and similarity to reference genomes. HPV was not identified as a major virus in our analysis of the retroauricular crease; however, molluscum contagiosum virus, a poxvirus that causes cutaneous growths that become severe in immunocompromised states, was present in high relative abundance in agreement with the original published findings (Oh et al., 2014).

Similar to our primary analysis, we identified 158 hypervariable loci within the Staphylococcus phage communities, and observed only 12 hypervariable loci associated with Propionibacterium phages, further highlighting the overall lack of genetic variability of the Propionibacterium phage communities. The Staphylococcus phage hypervariable loci were associated with purifying selection, yielding a pN/pS ratio slightly below 0.4 (Fig. 5B), recapitulating our findings in the primary dataset.

Staphylococcus phage nucleotide substitutions were associated with transitions between guanine and adenine residues, as we observed in our primary analysis (Fig. 5C). The ti/tv ratio of these loci was 1.02. The most common amino acid substitution was proline to leucine (Fig. 5D) and the substitution properties appeared loosely specific based on charge and polarity (Figs. 5E and 5F), reproducing our findings (Fig. 2).

We also evaluated the reproducibility of SNP position between identical but independently assembled genomic contigs of the two studies. We quantified the proportion of SNPs in our primary dataset that were also found at the same position in the secondary dataset. This revealed a median of approximately 50% overlap between datasets (Fig. 5G). As a control, we generated a simulated dataset using randomly assigned SNP positions instead of those determined experimentally. This yielded a significantly lower median of approximately 15% shared nucleotide SNP calls (Fig. 5G), suggesting that the observed SNP position is not random. These data indicate that our findings are consistent across different skin virome populations and techniques of collection and sequencing.