Characterization of prophages in bacterial genomes from the honey bee (Apis mellifera) gut microbiome

Bacterial genomes

Publicly available genomes (complete and fragmented) of bacteria originally isolated from the gut of the honey bee, A. mellifera, were downloaded from GenBank between March and July, 2020 (n = 181) (Table 1, Table S1). The nine core bacterial phylotypes are represented in this dataset (n = 151 total strains across 17 species), as were two bacterial pathogens (n = 30). Each genome was from a uniquely named strain, and in the case of duplicate strain names in NCBI, the most complete and recent genome was used.

Prophage identification

To detect putative prophages, we used a combination of VirSorter2 (v2.2.2; Guo et al., 2021, p. 2), CheckV (v.0.7.0 (Nayfach et al., 2021)) and VIBRANT (v1.2.1; Kieft, Zhou & Anantharaman, 2020). First, all genomes were analyzed with VirSorter2 (settings –include-groups “dsDNAphage,ssDNA,NCLDV,laviviridae”). Resulting viral regions were retained if they scored at least 0.5 for double stranded DNA phage (dsDNA phage, n = 539). Host genome regions flanking these viral sequences were then trimmed with the CheckV ’contamination’ command. To exclude highly degraded phages, trimmed sequences were retained only if they were at least 5 kb in length (n = 462). Of these resulting sequences, a region was considered a putative prophage if it scored at least 0.9 with VirSorter2 (n = 357). Additionally, those with VirSorter scores of 0.5–0.9 were further analyzed with VIBRANT and were retained as putative prophage if VIBRANT also classified these sequences as virus (n = 74). This resulted in a final set of 431 putative prophages.

Prophage abundance and prophage composition in honey bee bacterial isolates

Two metrics for prophage presence were assessed: the absolute abundance of prophages (total number of distinct regions) in a single bacterial genome and the portion of the bacterial genome that was composed of prophage (referred to as prophage composition). Differences in prophage abundance and composition among bacterial species were tested using Kruskal-Wallis tests. All graphs were visualized using ggplot2 and PNWColors in R (Wickham, 2016; Lawlor, 2020).

Larger host genomes often contain more prophage regions (Touchon, Bernheim & Rocha, 2016). We examined this pattern with host genome size vs. prophage abundance and host genome size vs. prophage composition using Spearman’s correlation with the R library package cocor (Diedenhofen & Musch, 2015) and the Meng, Rosenthal and Rubin’s z test (Meng, Rosenthal & Rubin, 1992; Myers & Sirois).

Functional annotation of prophage region genes and detection of intact phages

The 431 putative prophages were first annotated using prokka (v. 1.14.6) against the Prokaryotic Virus Remote Homologous Groups database (PHROGs) database (downloaded from http://millardlab.org/2021/11/21/phage-annotation-with-phrogs/ on July 20, 2022) (Seemann, 2014; Terzian et al., 2021). If an amino acid sequence had no hit to the PHROGs database (e value < 10⁻⁶), or was classified by PHROGs as “other” or “unknown function”, additional functions were then detected by aligning the amino acid sequences to hidden Markov profiles from the EggNOG 5.0 databases for bacteria, eukarya, archaea, and viruses (Huerta-Cepas et al., 2019) using hmmsearch (-E 0.001) with HMMER (v. 3.1b2) (hmmer.org, Eddy, 2011). Hits with bitscores above 30 were retained (Roux et al., 2016).

Functional composition of prophage genes were compared based on the COG category assignment within the EggNOG annotation files for bacteria, eukarya, archaea, and viruses databases. EggNOG hits were manually re-assigned to the category of “Phage-associated” if the descriptive EggNOG name included the words: phage, baseplate, capsid, integrase, tail, tape, lysozyme, portal, holin, N-Acetylmuramoyl-L-alanine amidase, transposase, virus, or viral. Sequences that did not match to anything in the EggNOG database, or that were classified as belonging to multiple COG categories, were re-classified as COG Category S, Function Unknown. Functional composition of prophages using the combined results of PHROGs and EggNOG were then compared among bacterial host species by summing the proteins of all prophages detected in genomes of that species and visualized with the packages ggplot2 and microshades (Wickham, 2016; Dahl et al., 2022).

To estimate intactness of prophage regions, trimmed sequences were predicted as intact using the following (i) those identified as intact by PHASTER (Arndt et al., 2016), (ii) those that encoded at least 3 or more different phage hallmark genes (see: cornerstone genes in Zhou et al., 2011) and additionally encoded an integrase and/or transposase.

Identification of unique putative prophages

To identify prophage regions likely belonging to the same population (boundary defined in Roux et al., 2019), genomes of the 431 putative prophages were dereplicated via pairwise average nucleotide identity (ANI) analysis using dRep (v3.0.0 (Olm et al., 2017)) (settings: –ignoreGenomeQuality -l 5000 -pa 0.95 -sa 0.95 -nc 0.85) with LASTn (v1080; Kiełbasa et al., 2011). Prophage pairs were considered to be from the same population if the ANI was above 95% over 85% of the shorter sequence (Roux et al., 2019). The representative prophage region for each population was selected with dRep. In one case, a population cluster contained two prophages that were identified from different bacterial host species. In that case, both prophage sequences were retained as representatives. As a result, a total of 237 putative prophage regions were unique to themselves and did not share high similarity with other sequences. The remaining 194 sequences had high similarity with at least one other prophage sequence, resulting in 66 distinct populations of phage. After dereplication, a total of 303 representative phage were used for downstream functional annotation and taxonomic analysis.

Viral clustering and classification

Phages lack a universally conserved, high-resolution phylogenetic marker gene to classify them. As such, prophages are often grouped into clusters via gene-sharing network-based approaches or whole genome sequence alignments to reference phage (Bolduc et al., 2017). Dereplicated prophages in this study (n = 303), along with 24,289 phage reference genomes downloaded from the INfrastructure for a PHAge REference Database (INPHARED) on January 25, 2023 (Cook et al., 2021) were first classified with gene-sharing network approaches using vConTACT2 (v0.9.22; –rel-mode Diamond –db ‘None’ –pcs-mode MCL –vcs-mode ClusterON). vContact2 (v0.9.22) groups phages into subclusters, indicating phage are likely to be within the same viral genera, and broader viral clusters that indicate subfamily relatedness (relatedness somewhere between family and genus level) (Bin Jang et al., 2019). The resulting networks were visualized using edge-weighted spring-embedded models in Cytoscape (Shannon et al., 2003; Bolduc et al., 2017; Bin Jang et al., 2019). Although vConTACT2 can provide taxonomic classification if target phage cluster with reference phage, for our dataset, few of the resulting clusters included references phages that enabled classification at any taxonomic level. Therefore, while vConTACT2 gene-sharing networks were used to assess the relationships among honey bee prophages, we ended up relying more on whole genome alignments to reference phage to assign broader taxonomy.

To assign taxonomy based on similarities to reference viral groups, we compared the average amino acid identity (AAI) of proteins between our dereplicated putative prophages and phage from a dereplicated Caudoviricetes database (amino acids downloaded from the INPHARED on November 5, 2022, taxonomic metadata downloaded from INPHARED on March 16, 2023) (Cook et al., 2021). The 16,253 reference genomes were first dereplicated with dRep (–ignoreGenomeQuality -p 32 -l 20000 -pa 0.90 –SkipSecondary inphared_derepout_90mash) resulting in 2,600 representative reference genomes. Amino acid sequences of proteins encoded by the reference were predicted using prodigal (default on each genome). A BLAST database was made of the reference proteins (makeblastdb). The honey bee gut prophages protein predictions (produced via prodigal, described above) were aligned to this reference database with BLASTp (E value 0.001). Taxonomic class was assigned at the most precise rank possible using each prophage’s AAI similarity to its top hitting reference phage. If a putative prophage and a reference phage shared greater than 70% AAI across 85% of proteins, it was classified as belonging to the same genus as the reference (Turner, Kropinski & Adriaenssens, 2021). If a phage did not meet this threshold but still shared at least 30% AAI over 50% of proteins, it was classified as the same family as the reference phage (Turner, Kropinski & Adriaenssens, 2021). If a prophage could not be classified using the above thresholds, but its top reference hit was to a Caudoviricetes class phage, it was classified broadly as Caudoviricetes. Any putative prophage with a top reference hit to a completely unclassified phage was labeled as Unclassified.

Distribution of prophages across core bacterial hosts

Bacterial host species significantly varied in both the number and composition of prophages across all bacterial isolates, both core and pathogens (Kruskal-Wallis, p < 0.0001; Fig. 1 & Table S2). A total of 269 high-confidence prophage regions were predicted across the 151 core bacterial isolates before dereplication. Bacterial genomes ranged in size from 1.3 Mb–3.5 Mb (Table S2). Among the core bacterial species, the number of prophages per isolate varied from 0 to 7 and prophage composition of the bacterial genome varied from 0–7.01% (Table S2).

The core bacteria with the highest median number of prophages and prophage composition was S. alvi (3.0 ± 1.46 prophages; 2.58% ± 1.40 prophage composition) and G. apicola (3.0 ± 1.58 prophages; 2.04% ± 1.48 prophage composition) (Figs. 1A and 1B, Table S2). Notably, while Bombella apis, Bartonella apis, and F. perrara had the second highest median number of prophages per genome, L. melliventris had the second highest prophage composition (Figs. 1A and 1B; Table S2). All other core species had a median number of 1 prophage or less and a prophage composition of less than 1% (Figs 1A and 1B, Table S2). These included Bifidobacterium asteroides, Lactobacillus mellifer, Lactobacillus apis, and Lactobacillus kimbladii, which had a median of 0 prophages per genome (Figs 1A and 1B; Table S2).

Of these 269 putative prophage regions, 90 were estimated to be intact. Bacterial hosts with the highest median number of estimated intact prophages were G. apicola (1.0 ± 1.1), S. alvi (1.0 ± 0.60) and L. melliventris (1 ± 0.58) (Fig. 1C; Table S2). All three isolates of L. melliventris, 71% of G. apicola isolates (n = 38) and 56% of S. alvi isolates (n = 34) possessed at least one estimated intact prophage (Table S3). The second highest median of intact prophage was found in L. kullabergensis (0.5 ± 0.71), with one of the two L. kullabergensis isolates possessing an intact prophage (Fig. 1C; Table S3). All other core species had a median of 0 intact prophages (Fig. 1C, Table S2). However, intact prophages did still occur in these isolates at lower frequencies (Table S3). Intact prophage composition varied slightly compared to prophage abundance, with L. melliventris having the highest composition (1.81% ± 1.25), followed by S. alvi (1.16% ± 1.09), L. kullabergensis (0.97% ± 1.37), and G. apicola (0.89% ± 1.31) (Fig. 1D; Table S2).

Distribution of prophages across pathogenic bacterial hosts

A total of 162 prophage regions were predicted across the 30 bacterial isolates from the two pathogens, P. larvae and M. plutonius, before dereplication. The bacterial isolate genomes ranged in size from 2.0 Mb–4.8 Mb (Table S2). The number of predicted prophage regions ranged from 1 to 22, and prophage composition ranged from 0.96%–13.85% of the bacterial genome (Table S2). P. larvae had the highest median number of prophages per genome (8.0 ± 5.33) and prophage composition (6.40% ± 3.08) across all bacterial species, including the core bacteria (Figs. 1A and 1B; Table S2). The median number of prophages found per genome of M. plutonius was 2.0 ± 0.64, and its median prophage composition was 2.77% ± 0.83 (Figs. 1A and 1B; Table S2).

Of the 13 P. larvae isolates analyzed, all had at least one intact prophage, with a median of 4 ± 1.98 intact prophage per isolate (TableS-Freq; Fig. 1C; Table S2). M. plutonius had a median of 1 ± 0.35 intact prophage, with 16 of 17 isolates possessing an intact prophage (Fig. 1C; Table S2; TableS-Freq). The total intact prophage composition of P. larvae (4.57% ± 3.10) was higher than M. plutonius (1.71% ± 0.97) (Fig. 1D; Table S2).

Prophage abundance is more correlated with bacterial genome size than prophage composition

Bacterial genome size was positively correlated with both the number of prophages (Spearman’s ρ = 0.55, p < = 1.6e−15, Fig. 2A) and prophage composition (Spearman’s ρ = 0.34, p = 2e−6, Fig. 2B). Meng, Rosenthal and Rubin’s correlation comparison determined that bacterial genome size had a stronger correlation with the number of prophages than with prophage composition (z = 5.9819, p-value < 0.001, 95% CI [0.18–0.35]) (Meng, Rosenthal & Rubin, 1992).

Prophage-encoded genes varied by bacterial host

Nearly half (48.4%) of all predicted protein-encoding genes of the dereplicated prophages were classified broadly as phage-associated. A total of 2.3% of all genes were classified by PHROGs specifically as “Moron, auxiliary metabolic genes, and host takeover” (Table S5). S. alvi had the highest frequency of genes in this category (66 genes, 3.3% frequency), due to its large collection of phage toxins and antitoxins (Fig. 3). These toxins shared homology to Doc-like toxins (10 genes), BstA abortive infection systems (five genes), HicB-like toxin-antitoxin systems (17 genes), and a RelE-like toxin (one gene), among other undescribed toxins (33 genes) (Table S7). P. larvae had the second highest number of genes (124 genes, 2.9% frequency), in the PHROGs category: moron, auxiliary metabolic genes, and host takeover (Fig. 3). P. larvae prophages possessed a variety of phage toxins, antitoxins, and additionally, genes associated with bacteriocin production or immunity (Table S7). Some P. larvae prophages also possessed a gene associated with quaternary ammonium compound-resistance proteins (two genes), or to ABC transporters (seven genes), indicating the presence of antimicrobial resistance genes (Table S7). One other antimicrobial-associated gene in a P. larvae prophage, originally classified as “unknown” by PHROGS, was annotated as a bacitracin resistance protein, BacA, by the eggNOG database (Table S7). Additionally, P. larvae prophages possessed genes for flavodoxin, a NifU-like Fe-S cluster assembly protein, and a phosphoadenosine phosphosulfate (PAPs) reductase (Table S7). The PAPs reductase also occurred frequently in prophages isolates from G. apicola (13 genes), and once in a B. asteroides prophage (Table S7).

The majority of remaining predicted protein-encoding genes (49.0%) were classified broadly as Function Unknown (COG Category S or no hits to the COG database) after failing to match to the PHROGs database (Table S5). Coding regions that were able to be classified broadly by COG were distributed into the categories of Information Storage and Processing (1.6%), Metabolism (0.5%), or Cellular Processing and Signaling (0.5%) (Table S4). Within these broader categories, COG category K (Transcription—Information Storage and Processing) occurred most frequently at 0.8% (Table S5). Amino Acid Transport and Metabolism (COG category E) and Carbohydrate Transport and Metabolism (COG category G) were the most commonly occurring metabolic categories, each corresponding to 0.1% of genes (Table S5, Table S7). Genes associated with amino acid metabolism were found in prophages of G. apicola (eight genes, 0.4% frequency) and P. larvae (five genes, 0.1% frequency) (Fig. 3, Table S6). Five different prophages of P. larvae appear to possess a gene homologous for a thermolysin metallopeptidase, while two share a gene associated with glycine amidinotransferase activity (Table S7). The amino acid genes associated with G. apicola prophages are more varied but appear to be associated with arginosuccinate and glutamate (Table S7). Prophage genes associated with carbohydrate metabolism were found in several bacterial hosts: M. plutonius (three genes, 0.8% frequency), G. apis (four genes, 0.7%), L. kunkeei (one gene, 0.41% frequency), Bartonella apis (one gene, 0.3% frequency), G. apicola (four genes, 0.2% frequency), and P. larvae (two genes, <0.1% frequency) (Fig. 3, Table S6). Interestingly, a single M. plutonius prophage possessed all three carbohydrate-associated proteins identified for that bacterial host: a mannitol dehydrogenase, a MFS/sugar transport protein, and a protein associated with the dehydration of D-mannotate (Table S7). Similarly, a single G. apicola prophage possessed all four of the carbohydrate-associated proteins identified for that bacterial host species: three proteins associated with the glycerate kinase family and one associated with alpha-ribazole phosphatase activity (Table S7). In contrast, three separate prophages of G. apis possessed a gene for alpha-galactosidase, with one also carrying a gene for an exopolysaccharide biosynthesis protein (Table S7). Additionally, two separate prophage species of P. larvae both possessed a gene associated with the D-gluconate metabolic process (Table S7). The remaining COG metabolic categories appear to be more rare, although they did still occur (Fig. 3, Table S6). For example, the sole B. indicum prophage that was identified possessed genes associated with nucleotide metabolism (a purine-cytosine permease) and secondary metabolite metabolism (isochorismatase) (Table S7).

Closely related prophages typically share host species

Using vConTACT2, a total of 208 honey bee prophages from the dereplicated prophages (n = 303) were grouped into 56 subclusters, with the remaining prophages detected as outliers (n = 82) or unable to be incorporated into the network (n = 6) (Fig. 4). The bacterial host strongly predicted clustering, with prophages of the same host species typically forming exclusive subfamily clusters (Fig. 4A). Only 10 of the 62 predicted subclusters included a prophage from more than one host, via either multiple honey bee bacterial species or the inclusion of non-honey bee reference phage (Fig. 4A). In several of these mixed clusters, the bacterial host genus was shared among isolates even if the species was not (Table S8).

In the resulting gene sharing network, honey bee prophages within the same clusters, and therefore often from the same bacterial host genome, were most closely positioned to each other (Fig. 4B). However, prophages from different host species sometimes appeared to share genetic similarities. For example, prophages from G. apicola and G. apis were closely positioned, as were some of the prophages from S. alvi and Bombella apis, and Bartonella apis (Fig. 4B). Interestingly, most of the prophages from Latobacillus spp., appear to share genes with each other, but prophages of L. kunkeei appear very diverse potentially due to few genes shared among these prophages (Fig. 4B). Furthermore, L. mellis prophages appear to share some similarity to prophages isolated from M. plutonius and P. larvae (Fig. 4B).

Most prophages from honey bee symbionts are unclassified Caudoviricetes

Of the 303 distinct prophage populations, 59.4% (180 prophage regions) were only able to be classified as Caudoviricetes, while an additional 31.7% prophages (96 prophage) were unable to be taxonomically classified using top reference hits (Fig. 5). Only the remaining 27 prophage regions (0.9%) were able to be classified at a family level or lower (Fig. 5). A total of 12 (16%) G. apicola prophages (n = 75) belonged to the family Peduoviridae (Class: Caudoviricetes), while one of three F. perrara prophages was predicted to belong to the Mesyanzhinovviridae family (Class: Caudoviricetes) (Fig. 5). The most common taxonomic ranks assigned to prophage from P. larvae isolates (n = 83) outside of Caudoviricetes (69 prophage, 66.3%) or unclassified (21 prophage, 20.2%) were the genera Fernvirus (seven prophage, 6.7%), Vegasvirus (three prophage, 2.9%) and Halyconevirus (two prophage, 1.9%) (Class: Caudoviricetes) (Fig. 5). A single P. larvae prophage each (1%) was classified as either the genera Lilyvirus or Dragolirvirus (Class: Caudovircetes).