Papillomavirus genome organization and conservation of the L1 jellyroll fold.

A. Prototypic arrangement of the papillomavirus (PV) genome. The PV genome consists of 6–8 open reading frames (green arrows) and an untranslated regulatory region (dark grey) arranged with conserved synteny. Core genes (dark sage), including the L1 major capsid protein (bright green) and its associated jellyroll region (light green), are found in all known PVs. The presence of accessory genes (light sage) varies between PV types. The PV capsid is composed of 72 self-assembling pentamers of the L1 protein, with the jellyroll region (represented here in both 2D topology and 3D structure) forming the principal structural fold. Protein and capsid structures retrieved from the HPV 16 L1 structure (PDB: 1DZL) and visualized in PyMOL37,38. B. The jellyroll fold is well-conserved across known PVs. A multiple sequence alignment (MUSCLE5 ‘super5’) from annotated PV L1 sequences in the Papillomavirus Episteme (PaVE, n = 614), clustered at 90% amino acid identity (n = 548, Supplementary File S1)27,78. Jellyroll beta-strands (B, C, D, E, F, G, H, I) were manually annotated based on the HPV 16 L1 structure (PDB: 1DZL) and residue positions in PyMOL37,38. Sequences are ordered by Jalview neighbour-joining tree (BLOSUM62) and visualized with the Taylor colour scheme106. Positional conservation is shown in the chartreuse histogram at the top of the alignment.

Discovery pipeline for novel papillomaviruses from Logan.

Sankey diagram illustrating the filtering of PV-homologous sequences retrieved from Logan. Cutoff values and filtering criteria are indicated in pink at each step. Open reading frames (ORFs) with homology to PV genes were identified with HMMs from PVDB2 using a bitscore threshold of >20, with L1 sequences further filtered to retain only those spanning the full jellyroll region. Graph traversal using PathRacer was used to retrieve additional full L1 contigs from libraries with partial hits (dashed teal line)31. To reduce redundancy and classify PVs at a ‘type’ level, full L1s were clustered at 90% nucleotide identity to generate type-representative sequences

Evolutionary, ecological, and geographic diversity of papillomaviruses.

A. Phylogenetic tree of PV L1 jellyroll protein sequences from Logan and NCBI Virus PV types. Sequences were aligned with MUSCLE5 ‘super5’ and the tree inferred with IQ-TREE2 using default parameters, the LG+G4 model, and 10,000 bootstrap iterations78,79. The tree was visualized with ggtree80. Grey branches represent known PV types; blue branches represent novel types identified in this study. Host organism groups (right) were manually assigned from SRA metadata, NCBI or PaVE annotations where available. Organism groups with fewer than nine representatives are classified as ’Other.’ B. Overlap between PV types identified from Logan and known PV types. Of 992 types in NCBI Virus, 646 (65%) were recovered by Logan. An additional 26 known types were identified through the NCBI nucleotide database (nt). A total of 383 novel types were identified. Note that for known NCBI Virus sequences, the known 646 are represented by 688 centroids from Logan. C. Library source of Logan PV sequences (unclustered). D. Distribution of known (grey) and novel (blue) PV types associated with host organism groups or E. ‘Other’ organisms. F. Geographic distribution of full-length L1+SRA libraries. Coordinates were inferred from BioSample metadata. Marginal histograms show the latitudinal and longitudinal density of detections.

Overlay of PV detections onto SRA sampling density.

Each grid square with any novel PVs detected are shaded teal to blue, while those with no novel types are shaded grey. Sampling effort is visualized from white to red, and is concentrated in the regions of North America, Western Europe, and East Asia. Statistically significant hotspots of novel type enrichment identified by Getis-Ord Gi* analysis (blue halos)34,35. Spatial autocorrelation of novel type proportion was confirmed by Moran’s I (I = 0.06, p < e⁻⁵)33.

Biome-level distribution of SRA sampling effort and PV discovery.

A. Geographic distribution of PV detections coloured by assigned biome. Each buffered grid square contains at least one library with a full-length L1 sequence, underlaid with the biomes of overall sampling efforts in the SRA. Biomes for each grid square were assigned based on the inferred geographic coordinates of the BioSample and geographical overlap with the WWF 2017 Ecoregions. The inset map shows the full global distribution of terrestrial and aquatic biomes used for classification. B. Total number of SRA sequencing runs per biome, ordered by sampling effort. Temperate Broadleaf and Mixed Forests account for the largest share of sequencing effort, followed by Ocean and Temperate Grasslands, Savannas and Shrublands. C. Total number of PV sequences detected per biome, ordered by sampling effort. Biomes are ordered identically in panels B and C to facilitate direct comparison. D. Proportion of novel PV types per biome, calculated as novel/(novel + known).

Case studies of novel papillomaviruses from Logan.

Two representative examples illustrating the biological depth recoverable from individual Logan contigs. A. White bellied pangolin-associated PV (contig SRR25256564_207500). B. Legless lizard-associated PV (contig SRR22028468_199653). For each case study, panels show (top) the genomic organization of the assembled contig with identified open reading frames, coloured as in Figure 1A; (middle) a phylogenetic subtree of closely related PV types with host associations indicated by colour, with grey numbers represent bootstrap support after 10,000 iterations, the inferred geographic origin of the source library (bottom left) and (bottom right) the predicted L1 protein structure generated by AlphaFold3, coloured by pLDDT confidence score88. The animal icons adjacent to each genome map correspond to the position markers on the full phylogenetic tree in Figure 3A.

Predicted gene content of five PV case studies.

For each PV, stop-stop ORFs were identified by NCBI ORFfinder, and by blastp against the nr database 32,107. For each gene, the top hit is shown with corresponding accession, e-value and amino acid percent identity (%ID). ORFs are visualized based on nucleotide coordinates of the assembled contig.

HMMs retrieved from Pfam for annotation of PV sequences from NCBI Virus.

Flowchart visualizing first round of PV searching, using PVDB1, and Logan v1.0.

Distribution of query coverage following blastn search of Logan L1 contigs.

Query coverage was calculated by dividing the matched residues over the total length of the input contig.

Proportion of PV gene coverage by SRA library type.

Stacked bars show the relative proportion of PV gene annotations (L1, L2, E1–E7), weighted by ka.f (k-mer-based coverage), across SRA library source categories grouped into DNA-, RNA-based, and miscellaneous. For multi-domain genes, ka.f values were averaged across E1 (E1_N, E1_DBD, E1_C) and E2 (E2_N, E2_C) sub-profiles. Chi-squared test was carried out between E and L ratios in pooled DNA and RNA sources.

Host group concordance of metadata-annotated PVs from Logan in comparison to NCBI annotations.

Each point represents one of the 646 NCBI Virus PV types found in Logan. The size and color represent the frequency at which the specific percent concordance and number of A. libraries or B. BioProjects were observed. The exact number and overall percent of types falling into each quadrant is indicated. Points with red diamond outlines indicate types in which more than 10 observations were made, but concordance was lower than expected. C, D. Distribution of associations derived from Logan metadata, where matches to NCBI Virus were found. On both a per-library and per-BioProject basis, the metadata-derived host for sequences with many observations but low concordance was visualized.

Host group enrichment analysis comparing Logan and NCBI Virus PVs.

Odds ratios (log scale) from Fisher’s exact tests comparing the proportion of each host group between PV types in Logan (n = 1,097) and NCBI Virus (n = 992). Odds ratios >1 (blue) indicate enrichment in Logan; odds ratios <1 (grey) indicate enrichment among NCBI references. Error bars represent 95% confidence intervals. Point size is proportional to the total number of PV types associated with each host group across both datasets. Significance levels are BH-adjusted for multiple comparisons (*p < 0.05, **p < 0.01, ns = not significant).

Green outlines of each PV-positive grid square represent those which have statistically significantly more novel PV enrichment following the 2000-iteration permutation test.

Visualization of all case studies.

A) Circular representation of the L1-based papillomavirus phylogenetic tree and generic papillomavirus genome. Blue tips of each branch indicate novel PVs uncovered in this study, and the colored outer ring corresponds to the associated host group. Generic PV genome was visualized using the expected average size of each ORF to show the consistent synteny throughout all case studies.

Best-hit Logan L1 nucleotide identities per host group.

Bars show the number of centroid clusters (n = 1,097), binned by 5% increments, and the percent nucleotide identity to the best hit in the nt database with at least 70% query coverage. Centroids with no hit are shown in the ‘no hit’ bin, and the dashed line marks the 90% species-level demarcation. The median bin for each host is colored in full saturation, and asterisks mark bins which are over-represented relative to overall frequency (Fisher’s exact test, Benjamini–Hochberg adjusted; * p < 0.05, ** p < 0.01, *** p < 0.001). Totals per host are annotated to the right.

Distribution of known and novel PV types across host organism groups.

Count of unique PV types associated with each host organism group, based on the sequences used for the phylogenetic analysis. Host groups were manually assigned from SRA metadata or PaVE annotation. Groups with fewer than nine total types are combined under the ‘Other’ category.