A diverse gut virome from Drosophila melanogaster

doi:10.21203/rs.3.rs-7226538/v1

Download PDF

Research Article

A diverse gut virome from Drosophila melanogaster

https://doi.org/10.21203/rs.3.rs-7226538/v1

This work is licensed under a CC BY 4.0 License

You are reading this latest preprint version

Background

Drosophila melanogaster is not only one of the most important models of antiviral immunity in invertebrates, but is also a powerful model for research of the gut microbiome. Although recent studies have continued to improve our knowledge of the fly gut microbiota, the viral component of the microbiome has remained unexplored.

Results

Here we explore the viral component of the Drosophila melanogaster gut microbiome using deep metagenomic DNA sequencing. We recovered 3035 phage sequences, resulting in 167 viral Metagenome-Assembled Genomes. The majority of these sequences are potentially novel bacteriophages from the order Caudovirales, which mainly target major gut bacteria of D. melanogaster, including Lactobacillus, Acetobacter, and Gluconobacter. Our functional annotation and discovery of auxiliary metabolic genes showed that these bacteriophages have the potential to influence microbial metabolism and genetic information processing. We also identified evidence of known fly pathogens Drosophila Kallithea nudivirus, Vesanto bidna-like virus, and Viltain densovirus, some of which were common in our studied populations.

Conclusions

Our findings reveal a complex and diverse phage community in the D. melanogaster gut microbiome, paving the way to study host-phage related research in the natural microbial communities.

gut microbiome

DNA virus

phage

Drosophila

The study of viruses, including bacteriophage (‘phages’; bacteria and archaea-infecting viruses) has recently garnered increasing interest in scientific communities due to the advancement of metagenomic techniques that facilitate the identification of viral genomes. Bacteriophages have been founded in nearly all investigated ecosystems and have been recognized as the most abundant and diverse biological entities on Earth [1]. They play an important role in microbial communities by affecting microbial diversity and metabolism [2, 3], and mediating phage-born antibiotic resistance within hosts [4].

Studies have also revealed the presence of auxiliary metabolic genes (AMGs) in bacteriophages [5]. These genes encode proteins that can regulate metabolic processes within host cells. In addition, bacteriophages have been found to regulate transcription and translation [6], and influence the decision to lyse or enter a lysogenic state [7], further underscoring the importance of characterizing bacteriophages and their interactions with their hosts.

In humans, it has been demonstrated that bacteriophage are mechanistically involved in human health and disease [8]. However, information about phages as an important component of microbiome infecting most species, is still lacking. This is also true for the well-studied model organism Drosophila, despite being a powerful model for gut microbiome research [9, 10] and the study of evolution and ecology of host-microbe interactions [11, 12].

The fruit fly D. melanogaster is a model organism that naturally interacts with a variety of microorganisms due to its feeding, mating, and ovipositing on fermenting and rotten fruits. Drosophila melanogaster has not only been established as a powerful model for gut microbiome research [9, 10], but is also well known as an important model for invertebrate antiviral immunity [13, 14]. However, most studies aimed at describing viruses in Drosophila have focused on the detection of RNA viruses of the host [15, 16], and studies on the diversity of DNA viruses have exclusively described eukaryotic DNA viruses [17]. Thus, almost nothing is known about the diversity or function of prokaryotic viruses—bacteriophages—associated with Drosophila. Here we aim to provide an initial characterization of the prokaryotic viral microbiome associated with D. melanogaster by applying metagenomic sequencing to screen the guts of wild-caught flies for viruses.

The metagenomic virome of D. melanogaster gut

We processed 851 metagenomes generated from D. melanogaster gut samples, resulting in nearly 38 billion high quality reads (range 1 to 90 million per sample). After mapping to the host, we retained 7.7 billion non-fly reads, and from these VIBRANT recovered 9124 metagenomic viral contigs (MVCs) that ranged in length from 1 to 545 kb with a mean length of approximately 7kb, and a total of 4964 contigs greater than 3Kb (Supplementary Table S1). After comparing the result with the finding from two other pieces of software (VirSorter2 and DeepVirfinder), filtering and dereplication, we identified 3040 unique contigs (Supplementary Table S2). Of the 3040 non-redundant contigs, 940 contigs were shorter than 3kb, and 1145 contigs were greater than 5kb. VIBRANT annotation identified 30 out of the 3040 nonredundant viral as complete circular genomes, 31 with high quality, 79 with medium-quality and 2900 as low-quality draft. Cenote-Taker 2 was also used, which checks for direct terminal repeats (DTRs) and inverted terminal repeats (ITRs), representing two common end features of complete viral genome. The Cenote-Taker 2 showed that 32 and 16 contigs out of 3040 contained DTRs and ITRs, respectively. These data showed 30 of 32 DTR-containing contigs matched with contigs identified as complete circular contigs by VIBRANT (Supplementary Table S2).

To further evaluate the final viral sequences, we applied checkV and accordingly most of the sequences (2821) were estimated to be low-quality genome fragments (< 50% complete), with only 67 with medium-quality (50–90% complete) and a total of 37 contigs were estimated to be high-quality (> 90% complete), of which 20 were predicted to be completed genome, also 115 contigs were not-determined (Supplementary Fig. S1). The chance of not-determined contigs to be false positive is low, as we retained them as viruses if they have been predicted with all three pieces of software (with VirSorter2 and DeepVirfinder quality scores ≥ 0.9). Therefore, the not-determined contigs are more likely to be too short to categorise, or the sequences are too divergent from the CheckV references—although Cenote-Taker 2 also could not predict halmarks for most of these 115 contigs. Next, we used VIBRANT to analyze whether the 3040 identified contigs were likely to have a lytic lifecycle. We found that 96.7% of them (2941 viral contigs) were predicted to have a lytic lifestyle, while only 3.3% represented a temperate lifestyle (Supplementary Table S2).

Protein Classification and Taxonomy of the virus contigs from the D. melanogaster gut

We evaluated the degree to which the newly recovered D. melanogaster gut-associated viruses were represented in available public databases (Fig. 1). Using a method that applies gene-sharing networks to classify sequences (vConTACT2), we determined that of 3040 non-redundant contigs over 1kb in length, 716 were singletons, 741 were outliers, 294 had genomes sharing overlap with multiple established clusters (i.e. cannot be unambiguously clustered), and 1289 were clustered (either with Refseq sequences or only with each other) (Supplementary Table S3 and Fig. S3A).

The 1289 clustered contigs formed 469 clusters, of which 434 clusters did not include any sequences from the Refseq database (Fig. 1A and Supplementary Table S3). The remaining 35 clusters, representing 67 contigs, were classified mainly as Siphoviridae, Myoviridae, and Podoviridae (Supplementary Fig. S2B). Furthermore, members of Schitoviridae, Herelleviridae, and Autographiviridae were present, although they comprised very small numbers of contigs (Fig. 1A and Supplementary Fig. S2B). In general, the resulting protein similarity network represented the newly identified phage sequences as broadly dispersed across the known viruses from the RefSeq database (Fig. 1A), despite the fact that 25% of the recovered sequences were also singletons (i.e. not presented in the network graph). The protein similarity network classification, however, showed that the higher-quality contigs (> 50% complete) tended to cluster more closely with previously identified viruses in the Refseq database, where almost 30% of these sequences assigned to known phages (Supplementary Fig. S3). Assignment success depended partly on contig length (Supplementary Fig. S4), with 53 (4.6%) of larger (> 5kb) contigs being assigned, but only 14 (~ 1%) of the smaller contigs (Supplementary Figs. S4A and S4B). However, the assignment of the large contigs to different phage families was comparable to the entire data set (Fig. 1 and Supplementary Figs. S2 versus Figs. S4).

To determine the family level taxonomy of the recovered contigs, we also characterised overlap in the composition of their viral protein families (VPF) with known viral taxa (Supplementary Table S4). Most of the recovered contigs have VPF compositions with family-level relatedness to Siphoviridae, Myoviridae, Podoviridae, and in lower ratios with Tectiviridae and Microviridae (an ssDNA virus). We also identified VPFs from archaeal virus families (Lipothrixviridae, Sphaerolipoviridae, Rudiviridae, and Turriviridae), NCLDV, and Nudiviridae—a family of DNA viruses known to infect eukaryotes, including Drosophila (Fig. 1B and Supplementary Table S4).

Of the 37 high-quality genomes, the VPF membership ratios of almost all suggest that they were related to bacteriophage families including Myoviridae, Siphoviridae, Podoviridae, and Microviridae (Fig. 1C). Although 27 (73%) of these almost complete phage genomes did not cluster with known reference phages (Supplementary Fig. S3), their VPF membership ratios were primarily related to Caudovirales (Fig. 1B and Supplementary Table S4), which may point to novel families or genera within this order.

Host Assignment of bacteriophage contigs

To predict the hosts of the recovered phage sequences, we took advantage of the VPFs as well as CRISPR-spacer sequences. Using the VPF approaches with a membership ratio and a confidence score cutoff of 0.6 and 0.8 respectively, 650 contigs could be linked to a host domain, of which almost all were assigned to the bacterial host (Supplementary Table S5). Further investigation showed that 267 of these phage contigs were linked to 15 host families and 154 of them could be related to 20 host genera. All identified potential hosts reflect either the gut or environmental bacteria (Supplementary Table S5). Of the 156 contigs that hit the host at genus-level, 43 could be related to the genus Pseudomonas, 30 to Lactococcus, 25 to Staphylococcus, 10 to Escherichia, and eight to Lactobacillus (Supplementary Table S5).

On the other hand, the CRISPR-spacer approach, which was used to predict bacterial hosts of determined phages, linked 488 bacteriophage contigs to 61 unique bacterial genera. Of these 488 contigs, most of them (> 300 contigs) hit against Drosophila gut specific bacteria (such as 74 Lactobacillus hits, 66 Gluconobacter hits, 48 Acetobacter hits, 31 Acinetobacter hits, 20 Pseudomonas hits, 19 Serratia hits, 15 Corynebacterium hits, 13 Klebsiella hits, 13 Komagataeibacter hits, etc.) (Supplementary Table S5).

Comparing the result from VPF and CRISPR-spacer approaches showed 29 contigs that hit the host at genus-level shared between two approaches, but only 69% of them (20 contigs) hit the same bacteria genus using both approaches (Supplementary Table S5). Further investigation showed that the VPF tools also predicted the presence of similar genus as CRISPR-spacer for these nine contigs, but were excluded due to our inclusion threshold (including only VPFs with a membership ratio and confident score of ≥ 0.6 and ≥ 0.8, respectively).

Functional Potential of D. melanogaster associated phage contigs from gut

Annotations from KEGG, Pfam, and VOG databases were used to gain insight into the functional potential encoded by the recovered contigs. Based on the available databases, only 50% (1549 contigs) of the contigs were assigned a function based on KEGG. However, the majority of phage associated functional genes were annotated with metabolism functions, including nucleotide, amino acid, lipid, energy and carbohydrate metabolism (Supplementary Table S6). Furthermore, functional genes associated with DNA repair, replication, and recombination, also transcription, were including quite large proportion (Supplementary Table S6).

We identified a total of 149 phage contigs that encoded at least one AMG and could be assigned to 42 unique KEGG orthologs (KO) (Supplementary Table S6). The most abundant AMGs (in terms of number of contigs involved) encoding dcm/DNMT1, folA/DHFR, mec, queE, NAMPT and folE/GCH1. These genes encode, respectively, DNA (cytosine-5)-methyltransferase 1, dihydrofolate reductase, the [CysO sulfur carrier protein] -S-L-cysteine hydrolase, the 7-cyano-7-deazaguanine synthase, nicotinamide phosphoribosyl transferase, and GTP cyclohydrolase IA (Supplementary Table 6). In other words, the AMGs derived from D. melanogaster tend to encode products mainly for amino acid, cofactor and vitamin metabolism and folding, sorting and degradation.

To further explore the functional potential of the identified contigs, we projected KO accession onto the KEGG pathways (Fig. 2). We used a total of 1549 phage contigs, which could be annotated based on KEGG database and assigned to 419 unique KO, for KEGG pathways projection (Supplementary Table 6). In some cases, the retrieved pathways represent functions that could directly affect bacteria, such as involvement in biofilm formation, quorum sensing, and the bacterial secretion system. In general, the retrieved pathways reflect a wide range of metabolic functions, including carbohydrate, energy, lipid, nucleotide, and amino acid metabolism. Furthermore, glycan biosynthesis and cofactors, vitamins, terpenoid and polyketide metabolism, and xenobiotic degradation were presented (Fig. 2). The identified genes in D. melanogaster gut phages are also involved in folding, sorting, degradation, transcription, and translation, as well as in replication and repair.

Genome binning of phage contigs from metagenomic data from D. melanogaster

We used vRhyme to construct viral metagenome-assembled genomes (vMAGs), resulting in the binning of 767 contigs into 188 vMAGs (Supplementary Table 7). However, after post-binning filtration based on potential contamination of vMAGs (bins that could include sequences from multiple viral genomes), in which we filtered out bins with ≥ 6 redundant proteins, we ended up with 167 vMAGs. We then checked the relative abundance of the final vMAG sets across the collected samples and observed variation related to the sampling locations (Fig. 3A). Out of 167 vMAGs, 11 were identified in more than 70% of the samples. Among these, two genomes, which their VPF membership ratios primarily belong to the Myoviridae family, were detected in all 851 samples (Fig. 3B and Supplementary Fig. S5).

Next, we further investigated the pattern of viral sequences across the samples, using principal coordinate analysis (PCoA), and clustering for sampling locations was assessed using Adonis2 and Betadisper tests. A significant effect from sampling location was observed not just based on the vMAG sequences but also based on unbinned viral sequences (2381 contigs) (Fig. 3C and Supplementary Fig. S6).

Drosophila melanogaster DNA viruses

During the analysis of virus protein composition, we identified a number of contigs deriving from nudiviruses, which are known to include multiple pathogens of Drosophila melanogaster (Wallace, Coffman et al., 2021). A blastn search against the NCBI viral genome database showed that these contigs were derived from Drosophila Kallithea nudivirus (NC_033829) with 99–100% coverage and a 99% identity. By mapping the high-quality short reads to the reference genome, we found that Kallithea virus occurs frequently, with a high overall prevalence of 54.2% (in 461 individuals of 851 samples with genome copy-number threshold ≥ 1% of the fly genome copy-number). The prevalence of Kallithea virus varied among four different sampling locations, ranging from 33.9% to 82.1%. We also detected Drosophila Esparto nudivirus, but its prevalence was extremely low, with only affected two flies.

By mapping reads to all previously-reported DNA virus genomes from Drosophila melanogaster [15, 17, 18], we identified two further previously-known pathogens, including Drosophila Linvill Road densovirus and Vesanto bidna-like virus. Drosophila Linvill Road densovirus was extremely rare, being detected in only seven individuals, all from sampling location 1, indicating a prevalence of 2 percent in that location and a total prevalence of 0.8% (95% confidence intervals 0.78%-0.82%). On the other hand, Drosophila Vesanto virus, a multi-segmented bidna-like virus, occurred at relatively high prevalence of 20% (95% confidence interval 18–23%; 172 samples out of 851) using a threshold of ≥ 1% (Fig. 4).

Almost all previous sequencing of Vesanto virus has been based on large metagenomic pools, leaving the associations between putative segments uncertain, as many independent infections are represented in each pool [17]. Here, based on single fly sequencing, we did not detect segments 9 or 12 —which were also absent from some population pools reported in [17]. Detection of the other segments was variable, with segment 5 being detected in 165 out of 851 samples (Fig. 4A, i.e. nearly all individuals infected with the Vesanto virus), segments 3 and 10 being detectable in 86 and 56 samples respectively and segments 1, 2, 4, 6, 7, and 8 each appearing in more than 35 individuals. Segment 11, was extremely rare, detected in only three samples.

Excluding segment 11, all other segments were detectable in those infected individuals that had the highest Vesanto virus read numbers (e.g., samples DE5797 and DE5682, which exhibited almost 2.5 and 1.8 million Vesanto virus reads, respectively). This is consistent with some segments (e.g. segments 9, 11, and 12) being ‘optional’ or ‘satellite’ components and other segments differing in copy-number, with segment 5 being the most readily detectable.

The virome is increasingly acknowledged as an essential part of the microbiome, playing a key role in coevolution with prokaryotic and eukaryotic hosts by both increasing and decreasing host fitness [19]. Here we used a metagenomic approach to elucidate the structure and function of the viral community associated with the gut of wild-caught Drosophila melanogaster, with a focus on prokaryotic viruses.

Phage associated with the Drosophila microbiome

Our study found that D. melanogaster harbors a diverse community of temperate and lytic bacteriophages in their gut (Supplementary Table S2). We also observed distinct variations in the structure of D. melanogaster bacteriophage communities from different sampling locations. This could suggest the gut microbiome in association with the host and environment is highly dynamic [20, 21].

We used two methods for host assignment of viral sequences: the VPF identification and the CRISPR spacers similarity method. Both methods assigned a relatively low percentage of the phages to a particular host (smaller than 16% in total). This is consistent with previous findings in honeybees [22] or humans [23]. Therefore, the number of viruses that could infect members of the fly core gut bacterial microbiome may be much higher than our results suggest. However, nearly 70% of the assigned phages were found to target specific Drosophila gut bacteria such as Lactobacillus, Acetobacter, and Gluconobacter [24–26]. This indicates that the phage community recovered is at least in part a true representative of the fly gut virome. Nevertheless, some retrieved sequences could represent food or components of the substrate-associated phage communities, as D. melanogaster feeds on microbes that grow on decaying substrate.

The majority of the viral sequences identified in this study potentially belonged to novel phage genera, with only a small percentage sharing amino acid identity with previously known phages. Almost 43% of identified sequences clustered into putative viral clusters, with only 5% deriving from a previously established reference genome based on ICTV (International Committee on Taxonomy of Viruses). The accuracy of the classification of the vConTACT2 algorithm [27] has been shown to be high, over 95% based on ICTV, supporting the classification of the sequences in this study. This finding indicates that the D. melanogaster microbiome is associated with a diverse range of phages, and that the current knowledge of phage diversity is limited. The fact that we lack full knowledge of phage diversity is further reinforced by similar observations in previous studies on honeybees [22] the human gut [28], and soil phages [29], which also found that only a small percentage of their sequences matched previously identified phages (less than 12%). These findings underscore the importance of continued efforts to study and understand the diversity of phages associated with different host organisms.

Phages are known to facilitate genetic material transfer in host communities, and can affect host phenotypes [30], suggesting that phage-associated proteins may play important roles not just in their life cycles but also in their host. The limited information available to categorize their functions highlights our incomplete understanding of these processes. Nevertheless, the prokaryotic viral-associated proteins that could be assigned to a function were predominantly related to metabolic and DNA processing pathways, in this study. This finding is consistent with observations in honeybees [22]. It indicates that the prokaryotic virome associated with the insect gut includes a vast range of functional genes related to metabolism, which may directly affect the gut ecosystem.

Phages are also shown to contribute to the evolution of their hosts [31], particularly in the rewiring of metabolic pathways to meet their requirements [32, 33]. Thus, the presence of important metabolic pathways (such as energy, vitamin, carbohydrate, and amino acid metabolisms) in the Drosophila-gut-associated virome, similar to observations in honeybees [22] and lizard [3], suggests the phages have the potential to modulate the metabolic state rather than solely rely on the host for resources. On the other hand, some pathways represented in the phages, particularly those related to genetic information processing and nucleotide metabolism, may reflect the phages' rewiring strategy.

In the context of phage-host interaction, the presence of accessory metabolic genes (AMGs) can determine the success of phage proliferation [5]. This is shown by the prevalence of the dcm gene, which encodes a methyltransferase that protects viral DNA from host restriction enzymes [34, 35]. Similarly, the 7-cyano-7-deazaguanine synthesis genes found also in the D. melanogaster bacteriophages offer phage protection against host defenses [36]. Our analysis shows that AMGs in D. melanogaster primarily encode products for amino acid, cofactor and vitamins and folding, sorting and degradation. These findings are consistent with previous studies on human and urban environments [2]. We speculate that these AMG functions not only reflect the needs of the phages, but also the environment and the host [29]. Such symbiotic interactions are critical to maintaining the balance of the gut ecosystem and sustaining the health of D. melanogaster individuals [37, 38].

Interestingly, our investigation also revealed the presence of pathways for biofilm formation, quorum sensing, and bacterial secretion system pathways within the D. melanogaster viral communities. These findings suggest that viruses may interfere with their bacterial hosts' processes and even influence other bacteria in the same ecosystem [22]. Moreover, the detection of secondary metabolites, terpenoids and polyketides metabolic pathways further indicates an intriguing interaction between phages and hosts. In light of these discoveries, we emphasize that bacteriophages play a critical role in microbial metabolism, and a balanced interaction between them is essential to sustain the continuum of the gut ecosystem. Ultimately, this intricate interplay between phages and bacteria indirectly influences the health and development of D. melanogaster individuals.

DNA viruses of Drosophila

DNA viruses of Drosophila melanogaster are thought to be relatively rare [17], and consistent with this we found only four viruses that are likely to infect the fly rather than its microbiome. These included Drosophila Kallithea nudivirus, Drosophila Esparto nudivirus, Drosophila Linvill Road densovirus, and Drosophila Vesanto virus.

Drosophila Kallithea nudivirus is a large dsDNA virus that is closely related to a number of other known Drosophila pathogens and has been shown to be pathogenic in adult flies [15, 17, 34, 39]. We found that Drosophila Kallithea nudivirus had high prevalence across our four sampling sites, with infection rates exceeding 50% of individuals. This is consistent a previous metagenomic study of DNA viruses infecting European D. melanogaster, which also reported that Drosophila Kallithea nudivirus infects more than half of populations [17].

Drosophila Vesanto virus was the second most common virus of Drosophila, being detected in more than 20% of the individuals. Here, using data from single flies, we detected only 10 out of 12 previously described segments of Drosophila Vesanto virus, and we observed distinct variation in the detection of the segments—with segment 11 appearing extremely rare. Based on large pool-sequencing datasets, [17] previously proposed three hypotheses to explain this pattern: that there are 12 segments with extreme variation in copy number; that there is a reassorting community of segments in which some are optional or satellites; or that Vesanto is not a single virus, and multiple viruses appeared in each of the pools presented there—perhaps infecting a component of the cellular microbiome. The very strong association we see between the majority of the segments (and the absence of shared eukaryotic microbiome from some affected individuals) strongly suggests that these are indeed segments of the same virus, as argued by [17] on the basis of a single infected lab line. The complete absence of two segments from all affected flies further confirms that some segments are indeed optional or satellites. However, the lack of negative correlations among most segments, along with the generally consistent pattern of copy-number variation, suggests that in general (apparent) absences are likely driven by the relatively low titre of some segments, rather than segments being optional or homologous being able to substitute for each other. Such variation in genome copy number among infections may hint that Vesanto virus has a multipartite lifestyle [40, 41].

Sample Collection and Metagenomic Sequencing

A total of 851 male Drosophila melanogaster were collected from four sampling locations in Germany (Baden-Württemberg: for more details see Supplementary Fig. S7), over a two-month period between late August and early October-2019. For sample collection, we used traps baited with a mixture of apple and cherry, following a standardized protocol [18]. The flies were anesthetized on a fly pad using CO₂ and subsequently transferred into empty Eppendorf tubes. The Eppendorf tubes were then dipped into liquid nitrogen to freeze the samples. Following this, all fly samples were dissected to extract the gut, and then the individual gut samples were stored in PBS at -80°C until DNA extraction.

We extracted DNA from each gut sample individually following the DNA extraction protocol described in detail in [18] with slight modification. Briefly, the samples were homogenized using a bead beater (Qiagen Tissue Lyzer II), with the addition of lysozyme to improve the efficiency of nucleic acid extraction. Protein was digested with Proteinase K, and RNA was removed using RNAse. We precipitated DNA using a combination of Phenol-Chloroform-Isoamyl alcohol and after washing with alcohol the DNA was resuspended in nuclease-free water. The NEBNext® Multiplex Oligos with unique dual index primer pairs for Illumina was used for library preparation. Then samples were sequenced on the Novaseq 6000 platform with 150bp paired end sequencing strategy at the Life & Brain research centre (University Hospital Bonn).

Metagenome assembly

We used Cutadapt [42] and FastQC [43] implemented in TrimGalore v0.6.6 [42] to trim adapters and control the quality of the raw reads. Then, to remove human and host sequences, we mapped the genomic DNA of the data from each of the 851 sample to the host D. melanogaster and the human genome (hg38), using bbmap [44]. Finally, we assembled the remaining paired-end reads using MEGAHIT [45]. To improve the assembly of low-abundance components of the microbiome, we co-assembled reads from multiple individuals, pooling by sampling location, for a total of four co-assemblies.

Recovering high-quality viruses and Binning of viral genomes

We predicted viral contigs among metagenome contigs using VIBRANT v1.2.1 [46] with the default parameters. VIBRANT uses machine learning and a protein similarity approach to accurately recover bacteriophages [47]. In outline, VIBRANT uses the contigs predicted by Prodigal v2.6.3 [48] and annotates them with the Kyoto Encyclopedia of Genes and Genomes (KEGG) KoFam [49], Pfam [50], and virus orthologous group (VOG) [51] databases using HMMER v3.3.2 [52]. At the same time, we also inferred auxiliary metabolic genes (AMG) using VIBRANT. VIBRANT employs neural networks of protein signatures, and a novel v-score metric to overcome conventional limitations, enabling comprehensive identification of lytic viral genomes, thereby capable of providing estimates of genome quality and differentiation between lytic and lysogenic viruses.

The completeness and quality of all phage contigs were assessed using CheckV v1.4 [53] with the parameters “end_to_end” and removing all contigs shorter than 3kb. To ensure the quality of retained phage contigs, we ran VirSorter v2.2.3 [54] and DeepVirfinder (v. 1.0) [55] against the same assembled data. We considered the phage contigs recovered by VIBRANT as true bacteriophages if the prediction scores of VirSorter2 (groups include-dsDNAphage, ssDNA viruses, NCLDV-Nucleocytoviricota) and DeepVirfinder also were above ≥ 0.5.

DeepVirFinder has especially been advocated for the identification of short (1 < kbp) phage genome fragments [56]. The dsDNA, ssDNA and NCLDV, respectively, represent double-strand and single-strand DNA viruses and nucleo-cytoplasmic large DNA viruses. However, to reduce the chance of false positives given the sensitive cut-off point for the length of the VIBRANT and VirSorter2 fragment of 3 kbp, the short contigs < 3 kbp) were also retained if they had been recovered by all three pieces of software and had VirSorter2 and DeepVirfinder quality scores ≥ 0.9. Contigs with provirus contamination were also retained if they exhibited a clear shift from the host genome and met the aforementioned criteria, after extraction of the proviral region.

All identified contigs from the four co-assemblies were pooled and dereplicated (–method longest) using vRhyme v.1.1.0 [57]. The quality of the final nonredundant phage sequences has been assessed using CheckV v1.4 [53]. We also used Cenote-Taker 2 v2.1.5 with “-am True” [8] to further assess and annotate the nonredundant recovered sequences.

Finally, we used the nonredundant viral contigs binning strategy to group scaffolds into putative genomes using vRhyme [57] with the default parameters. As part of bin construction, vRhyme checks for protein redundancy and this metric can be used for post-binning filtration (Redundant protein < 1 unlikely to be contaminated: 2 < 5 may not be contaminated; >6 are more likely to be contaminated) [57].

To perform a quantitative analysis of the obtained phage contigs, we used metaWRAP with “quant_bins” module to map the high-quality reads against the recovered contigs. Then the average abundance of each (viral metagenome-assembled genome) vMAG in each of 851 sample was calculated using Salmon [58], reporting quantities as transcripts per million (TPM) [59] and taking the length-weighted average of the contig abundances.

We performed principal coordinate analysis (PCoA) and visualized the Bray–Curtis distances of relative viral abundances at each sampling location using the vegan package (v2.5-7) in R v4.2.0 [60]. Then the observed variation was further tested using Adonis2, also Betadisper test in the vegan package [61], to gain a better understanding of the underlying patterns and structure in the data set as well as validate the findings obtained through PCoA. Adonis2 helps determine if the observed groupings in PCoA are statistically significant, while Betadisper provides insights into the within-group variability and helps identify if the groups identified by PCoA have significantly different spread or variability. We also made a heatmap to visualize the differences between phage composition of different localities using Seaborn clustermap v 0.12.2 [62].

Protein clustering and taxonomy of viral contigs

We evaluated viral contig similarity and constructed a viral protein clustering network with vConTACT2 0.11.3 [27] using the viral Refseq v211 database. The vConTACT2 method uses gene-sharing networks to assign the taxonomy of viruses based on their sequence. The viral sequences are then classified as ‘clustered’ (high-confidence clustering), ‘overlap’ (sharing overlap gene content), ‘outlier ‘(weakly associated with a cluster) or ‘singleton’ (few or no shared gene content) [27]. The network was visualised using Cytoscape v3.9.1 [63]. We performed host assignment and taxonomic annotation at the family and genus level with VPF-Class [64] using VPF classified files from 17 August 2022 to build the reference database. For preliminary taxonomic annotation, we considered a VPF membership ratio cutoff of 0.2 and 0.3 at the family and genus level, respectively.

Host Assignment of recovered bacteriophage contigs

To assign viral genomic sequences to their specific host, we applied two strategies. First, we used VPF-class tool, which takes advantage of the Viral Protein Families (VPFs) and creates a protein-based database from the IMG/VR database. The VPF-class tools then conducts a protein family search for each query genome and concludes the host prediction based on the distribution of these VPFs in reference phage genomes, first at the domain and then at the host family and genus levels [64]. In the second approach, we used CRISPR spacer-based bacterial host predictions strategy [65]. The CRISPR spacer-based approach uses information from CRISPR–Cas, which works as adaptive immune systems, and has access to a database of more than 11 million spacers [65]. We therefore searched the spacer database with our viral sequences using blastn and predicted bacterial hosts at the genus level on the basis of sequence similarity.

DNA viruses of Drosophila melanogaster

Initial taxonomic assignment of the recovered contigs suggested that at least one nudivirus was present. We therefore used blastn to test whether this was a previously known viral pathogen of Drosophila. To detect the presence of other previously-known DNA viruses of Drosophila, we mapped the high-quality short reads to the reference genomes of Drosophila Kallithea nudivirus, Drosophila Viltain densovirus, Drosophila Linvill Road densovirus, and Drosophila Vesanto bidna-like virus to quantify the amount of each virus in each fly [15, 17, 18]. We set a threshold of 1 per cent of the fly genome copy number to define infection status, and using these data, we were able to estimate the prevalence of these fly pathogens in our 851 samples.

AMGs

auxiliary metabolic genes MVCs:metagenomic viral contigs

vMAGs

viral metagenome-assembled genomes VPF:viral protein families

DTRs

direct terminal repeats ITRs:inverted terminal repeats

PCoA

principal coordinate analysis

Ethics approval and consent to participate

Not Applicable

Consent for publication

Not Applicable

Competing interests

All authors declare no conflicts of interest and have approved the manuscript for publication.

Funding

This work was supported by the German Research Foundation (DFG) with– Project number 408908608.

Author Contribution

FS and NA collected samples, and NA performed dissections. MHA carried out the laboratory protocols and conducted the bioinformatic analyses. MHA and DJO performed the statistical analyses. MHA drafted the manuscript with substantial contributions from DJO and input from all co-authors. All authors read and approved the final version of the manuscript.

Acknowledgement

The authors gratefully acknowledge the support of the European Drosophila Population Genomics Consortium (DrosEU), funded by a Special Topics Network (STN) grant from the European Society of Evolutionary Biology (ESEB). Additionally, we appreciate the support from the German Research Foundation (DFG) through project funding. We also thank the state of Baden-Württemberg for their support through bwHPC. The authors extend their gratitude to Dr. Mehregan Ebrahimi for his assistance in data analysis. Special thanks are also given to the anonymous reviewers for their valuable comments on the manuscript.

Data Availability

The sequencing data generated and supporting the findings of this study are deposited in the NCBI Sequence Read Archive (SRA) under BioProject accession number PRJNA1097789[ [https://www.ncbi.nlm.nih.gov/bioproject/PRJNA1097789](https:/www.ncbi.nlm.nih.gov/bioproject/PRJNA1097789) ]. Scripts for the bioinformatics processing are publicly available at [https://github.com/Minamehr/Manuscripts/tree/main/Virus-Drosophila-melanogaster](https:/github.com/Minamehr/Manuscripts/tree/main/Virus-Drosophila-melanogaster) .

Dion MB, Oechslin F, Moineau S. 2020 Phage diversity, genomics and phylogeny. Nature Reviews Microbiology. 18, 125–138. (10.1038/s41579-019-0311-5).
Kieft K, Zhou Z, Anderson RE, Buchan A, Campbell BJ, Hallam SJ, Hess M, Sullivan MB, Walsh DA, Roux S, Anantharaman K. 2021 Ecology of inorganic sulfur auxiliary metabolism in widespread bacteriophages. Nature Communications. 12, 3503. (10.1038/s41467-021-23698-5).
Luo X-Q, Wang P, Li J-L, Ahmad M, Duan L, Yin L-Z, Deng Q-Q, Fang B-Z, Li S-H, Li W-J. 2022 Viral community-wide auxiliary metabolic genes differ by lifestyles, habitats, and hosts. Microbiome. 10, 190. (10.1186/s40168-022-01384-y).
Moon K, Jeon JH, Kang I, Park KS, Lee K, Cha C-J, Lee SH, Cho J-C. 2020 Freshwater viral metagenome reveals novel and functional phage-borne antibiotic resistance genes. Microbiome. 8, 75. (10.1186/s40168-020-00863-4).
Thompson LR, Zeng Q, Kelly L, Huang KH, Singer AU, Stubbe J, Chisholm SW. 2011 Phage auxiliary metabolic genes and the redirection of cyanobacterial host carbon metabolism. Proceedings of the National Academy of Sciences. 108, E757-E764. (10.1073/pnas.1102164108).
Liu X, Hong T, Parameswaran S, Ernst K, Marazzi I, Weirauch MT, Bass F. J. I. 2020 Human Virus Transcriptional Regulators. Cell. 182, 24–37. (10.1016/j.cell.2020.06.023).
Dou C, Xiong J, Gu Y, Yin K, Wang J, Hu Y, Zhou D, Fu X, Qi S, Zhu X et al. 2018 Structural and functional insights into the regulation of the lysis–lysogeny decision in viral communities. Nature Microbiology. 3, 1285–1294. (10.1038/s41564-018-0259-7).
Tisza MJ, Buck CB. 2021 A catalog of tens of thousands of viruses from human metagenomes reveals hidden associations with chronic diseases. Proceedings of the National Academy of Sciences. 118, e2023202118. (10.1073/pnas.2023202118).
Ludington WB, Ja WW. 2020 Drosophila as a model for the gut microbiome. PLOS Pathogens. 16, e1008398. (10.1371/journal.ppat.1008398).
Younes S, Al-Sulaiti A, Nasser EAA, Najjar H, Kamareddine L. 2020 Drosophila as a Model Organism in Host–Pathogen Interaction Studies. Frontiers in Cellular and Infection Microbiology. 10, (10.3389/fcimb.2020.00214).
Ahlawat N, Geeta Arun M, Maggu K, Jigisha, Singh A, Prasad NG. 2022 Drosophila melanogaster hosts coevolving with Pseudomonas entomophila pathogen show sex-specific patterns of local adaptation. BMC Ecology and Evolution. 22, 77. (10.1186/s12862-022-02031-8).
Henry LP, Ayroles JF. Drosophila melanogaster microbiome is shaped by strict filtering and neutrality along a latitudinal cline. Mol Ecol. 2022;31:5861–71. https://doi.org/10.1111/mec.16692).
Xu J, Cherry S. Viruses and antiviral immunity in Drosophila. Dev Comp Immunol. 2014;42:67–84. https://doi.org/10.1016/j.dci.2013.05.002).
Mussabekova A, Daeffler L, Imler J-L. 2017 Innate and intrinsic antiviral immunity in Drosophila. Cellular and Molecular Life Sciences. 74, 2039–2054. (10.1007/s00018-017-2453-9).
Webster CL, Waldron FM, Robertson S, Crowson D, Ferrari G, Quintana JF, Brouqui J-M, Bayne EH, Longdon B, Buck AH. 2015 The discovery, distribution, and evolution of viruses associated with Drosophila melanogaster. PLoS Biol. 13, e1002210.
Webster CL, Longdon B, Lewis SH, Obbard DJ. 2016 Twenty-Five New Viruses Associated with the Drosophilidae (Diptera). Evolutionary Bioinformatics. 12s2, EBO.S39454. (10.4137/ebo.S39454).
Wallace MA, Coffman KA, Gilbert C, Ravindran S, Albery GF, Abbott J, Argyridou E, Bellosta P, Betancourt AJ, Colinet H. The discovery, distribution, and diversity of DNA viruses associated with Drosophila melanogaster in Europe. Virus Evol. 2021;7:veab031.
Kapun M, Barrón MG, Staubach F, Obbard DJ, Wiberg RAW, Vieira J, Goubert C, Rota-Stabelli O, Kankare M, Bogaerts-Márquez M et al. 2020 Genomic Analysis of European Drosophila melanogaster Populations Reveals Longitudinal Structure, Continent-Wide Selection, and Previously Unknown DNA Viruses. Molecular Biology and Evolution. 37, 2661–2678. (10.1093/molbev/msaa120).
Cadwell K. 2015 The virome in host health and disease. Immunity. 42, 805–813. (10.1016/j.immuni.2015.05.003).
Rothschild D, Weissbrod O, Barkan E, Kurilshikov A, Korem T, Zeevi D, Costea PI, Godneva A, Kalka IN, Bar N. 2018 Environment dominates over host genetics in shaping human gut microbiota. Nature 555, 210–5.
Gogarten JF, Rühlemann M, Archie E, Tung J, Akoua-Koffi C, Bang C, Deschner T, Muyembe-Tamfun J-J, Robbins MM, Schubert G. 2021 Primate phageomes are structured by superhost phylogeny and environment. Proceedings of the National Academy of Sciences. 118, e2013535118.
Deboutte W, Beller L, Yinda CK, Maes P, de Graaf DC, Matthijnssens J. 2020 Honey-bee–associated prokaryotic viral communities reveal wide viral diversity and a profound metabolic coding potential. Proceedings of the National Academy of Sciences. 117, 10511–10519. (10.1073/pnas.1921859117).
Shkoporov AN, Clooney AG, Sutton TD, Ryan FJ, Daly KM, Nolan JA, McDonnell SA, Khokhlova EV, Draper LA, Forde A. The human gut virome is highly diverse, stable, and individual specific. Cell Host Microbe. 2019;26:527–41. e525.
Staubach F, Baines JF, Künzel S, Bik EM, Petrov DA. 2013 Host Species and Environmental Effects on Bacterial Communities Associated with Drosophila in the Laboratory and in the Natural Environment. PLOS ONE. 8, e70749. (10.1371/journal.pone.0070749).
Broderick NA, Lemaitre B. Gut-associated microbes of Drosophila melanogaster. Gut Microbes. 2012;3:307–21.
Pais IS, Valente RS, Sporniak M, Teixeira L. 2018 Drosophila melanogaster establishes a species-specific mutualistic interaction with stable gut-colonizing bacteria. PLoS Biol. 16, e2005710.
Bin Jang H, Bolduc B, Zablocki O, Kuhn JH, Roux S, Adriaenssens EM, Brister JR, Kropinski AM, Krupovic M, Lavigne R. 2019 Taxonomic assignment of uncultivated prokaryotic virus genomes is enabled by gene-sharing networks. Nat Biotechnol. 37, 632–9.
Gregory AC, Zablocki O, Zayed AA, Howell A, Bolduc B, Sullivan MB. The gut virome database reveals age-dependent patterns of virome diversity in the human gut. Cell Host Microbe. 2020;28:724–40. e728.
Busse L, Tisza M, DiRuggiero J. 2022 Viruses Ubiquity and Diversity in Atacama Desert Endolithic Communities. Viruses. 14, 1983.
Virgin HW. The virome in mammalian physiology and disease. Cell. 2014;157:142–50.
Boyd EF. Bacteriophage-encoded bacterial virulence factors and phage–pathogenicity island interactions. Adv Virus Res. 2012;82:91–118.
Enav H, Mandel-Gutfreund Y, Béjà O. Comparative metagenomic analyses reveal viral-induced shifts of host metabolism towards nucleotide biosynthesis. Microbiome. 2014;2:1–12.
De Smet J, Zimmermann M, Kogadeeva M, Ceyssens P-J, Vermaelen W, Blasdel B, Jang B, Sauer H, Lavigne U, R. High coverage metabolomics analysis reveals phage-specific alterations to Pseudomonas aeruginosa physiology during infection. ISME J. 2016;10:1823–35.
Palmer WH, Medd NC, Beard PM, Obbard DJ. 2018 Isolation of a natural DNA virus of Drosophila melanogaster, and characterisation of host resistance and immune responses. PLoS Pathog. 14, e1007050.
Boye E, Løbner-Olesen A. The role of dam methyltransferase in the control of DNA replication in E. coli. Cell. 1990;62:981–9.
Hutinet G, Kot W, Cui L, Hillebrand R, Balamkundu S, Gnanakalai S, Neelakandan R, Carstens AB, Lui F, Tremblay C. D. 2019 7-Deazaguanine modifications protect phage DNA from host restriction systems. Nat Commun. 10, 5442.
Fernández L, Rodríguez A, García P. Phage or foe: an insight into the impact of viral predation on microbial communities. ISME J. 2018;12:1171–9.
Hwang Y, Rahlff J, Schulze-Makuch D, Schloter M, Probst AJ. 2021 Diverse Viruses Carrying Genes for Microbial Extremotolerance in the Atacama Desert Hyperarid Soil. mSystems. 6, e00385-00321. (10.1128/mSystems.00385-21).
Unckless RL. A DNA virus of Drosophila. PLoS ONE. 2011;6:e26564.
Sicard A, Yvon M, Timchenko T, Gronenborn B, Michalakis Y, Gutierrez S, Blanc S. 2013 Gene copy number is differentially regulated in a multipartite virus. Nature Communications. 4, 2248. (10.1038/ncomms3248).
Michalakis Y, Blanc S. 2020 The Curious Strategy of Multipartite Viruses. Annual Review of Virology. 7, 203–218. (10.1146/annurev-virology-010220-063346).
Martin M. 2011 Cutadapt removes adapter sequences from high-throughput sequencing reads. 2011. 17, 3. (10.14806/ej.17.1.200).
Wingett SW, Andrews S. 2018 FastQ Screen: A tool for multi-genome mapping and quality control. F1000Res. 7, 1338. (10.12688/f1000research.15931.2).
Bushnell B. BBMap: a fast, accurate, splice-aware aligner. Berkeley, CA: Lawrence Berkeley National Lab.(LBNL); 2014. (United States).
Li D, Luo R, Liu C-M, Leung C-M, Ting H-F, Sadakane K, Yamashita H, Lam T-W. 2016 MEGAHIT v1.0: A fast and scalable metagenome assembler driven by advanced methodologies and community practices. Methods. 102, 3–11. (https://doi.org/10.1016/j.ymeth.2016.02.020).
Kieft K, Zhou Z, Anantharaman K. 2020 VIBRANT: automated recovery, annotation and curation of microbial viruses, and evaluation of viral community function from genomic sequences. Microbiome. 8, 90. (10.1186/s40168-020-00867-0).
Heyerhoff B, Engelen B, Bunse C. 2022 Auxiliary Metabolic Gene Functions in Pelagic and Benthic Viruses of the Baltic Sea. Front Microbiol. 13, 863620. (10.3389/fmicb.2022.863620).
Hyatt D, Chen G-L, LoCascio PF, Land ML, Larimer FW, Hauser LJ. 2010 Prodigal: prokaryotic gene recognition and translation initiation site identification. BMC Bioinformatics. 11, 119. (10.1186/1471-2105-11-119).
Aramaki T, Blanc-Mathieu R, Endo H, Ohkubo K, Kanehisa M, Goto S, Ogata H. 2019 KofamKOALA: KEGG Ortholog assignment based on profile HMM and adaptive score threshold. Bioinformatics. 36, 2251–2252. (10.1093/bioinformatics/btz859).
El-Gebali S, Mistry J, Bateman A, Eddy SR, Luciani A, Potter SC, Qureshi M, Richardson LJ, Salazar GA, Smart A et al. 2018 The Pfam protein families database in 2019. Nucleic Acids Research. 47, D427-D432. (10.1093/nar/gky995).
Grazziotin AL, Koonin EV, Kristensen DM. 2016 Prokaryotic Virus Orthologous Groups (pVOGs): a resource for comparative genomics and protein family annotation. Nucleic Acids Research. 45, D491-D498. (10.1093/nar/gkw975).
Eddy SR. 2011 Accelerated Profile HMM Searches. PLOS Computational Biology. 7, e1002195. (10.1371/journal.pcbi.1002195).
Nayfach S, Camargo AP, Schulz F, Eloe-Fadrosh E, Roux S, Kyrpides NC. 2021 CheckV assesses the quality and completeness of metagenome-assembled viral genomes. Nature Biotechnology. 39, 578–585. (10.1038/s41587-020-00774-7).
Guo J, Bolduc B, Zayed AA, Varsani A, Dominguez-Huerta G, Delmont TO, Pratama AA, Gazitúa MC, Vik D, Sullivan MB, Roux S. 2021 VirSorter2: a multi-classifier, expert-guided approach to detect diverse DNA and RNA viruses. Microbiome. 9, 37. (10.1186/s40168-020-00990-y).
Ren J, Song K, Deng C, Ahlgren NA, Fuhrman JA, Li Y, Xie X, Poplin R, Sun F. 2020 Identifying viruses from metagenomic data using deep learning. Quantitative Biology. 8, 64–77. (10.1007/s40484-019-0187-4).
Pratama AA, Bolduc B, Zayed AA, Zhong ZP, Guo J, Vik DR, Gazitúa MC, Wainaina JM, Roux S, Sullivan MB. 2021 Expanding standards in viromics: in silico evaluation of dsDNA viral genome identification, classification, and auxiliary metabolic gene curation. PeerJ. 9, e11447. (10.7717/peerj.11447).
Kieft K, Adams A, Salamzade R, Kalan L, Anantharaman K. 2022 vRhyme enables binning of viral genomes from metagenomes. Nucleic Acids Research. 50, e83-e83. (10.1093/nar/gkac341).
Patro R, Duggal G, Love MI, Irizarry RA, Kingsford C. 2017 Salmon provides fast and bias-aware quantification of transcript expression. Nature Methods. 14, 417–419. (10.1038/nmeth.4197).
Wagner GP, Kin K, Lynch VJ. 2012 Measurement of mRNA abundance using RNA-seq data: RPKM measure is inconsistent among samples. Theory in Biosciences. 131, 281–285. (10.1007/s12064-012-0162-3).
R Core Team. 2014 R: A language and environment for statistical computing. MSOR connections. 1.
Oksanen J, Blanchet FG, Friendly M, Kindt R, Legendre P, McGlinn D, Minchin PR, O'hara R, Simpson GL, Solymos P. 2019 Package ‘vegan’. Community ecology package, version. 2.
Waskom ML. Seaborn: statistical data visualization. J Open Source Softw. 2021;6:3021.
Saito R, Smoot ME, Ono K, Ruscheinski J, Wang P-L, Lotia S, Pico AR, Bader GD, Ideker T. 2012 A travel guide to Cytoscape plugins. Nat Methods. 9, 1069–76.
Pons JC, Paez-Espino D, Riera G, Ivanova N, Kyrpides NC, Llabrés M. 2021 VPF-Class: taxonomic assignment and host prediction of uncultivated viruses based on viral protein families. Bioinformatics 37, 1805–13.
Dion MB, Plante P-L, Zufferey E, Shah SA, Corbeil J, Moineau S. 2021 Streamlining CRISPR spacer-based bacterial host predictions to decipher the viral dark matter. Nucleic Acids Res. 49, 3127–38.

No competing interests reported.

Download PDF

Editorial decision: Revision requested
25 Dec, 2025
Reviews received at journal
25 Dec, 2025
Reviews received at journal
22 Dec, 2025
Reviewers agreed at journal
14 Dec, 2025
Reviewers agreed at journal
13 Dec, 2025
Reviews received at journal
26 Oct, 2025
Reviewers agreed at journal
22 Sep, 2025
Reviewers agreed at journal
20 Sep, 2025
Reviewers invited by journal
18 Sep, 2025
Editor assigned by journal
11 Sep, 2025
Submission checks completed at journal
11 Sep, 2025
First submitted to journal
07 Sep, 2025

You are reading this latest preprint version

A diverse gut virome from Drosophila melanogaster

Status:

Version 1

Abstract

Background

Results

Conclusions

Figures

Background

Results

Host Assignment of bacteriophage contigs

Discussion

Methods

Sample Collection and Metagenomic Sequencing

Metagenome assembly

Recovering high-quality viruses and Binning of viral genomes

Protein clustering and taxonomy of viral contigs

Host Assignment of recovered bacteriophage contigs

Abbreviations

Declarations

Ethics approval and consent to participate

Consent for publication

Funding

Author Contribution

Acknowledgement

Data Availability

References

Additional Declarations

Supplementary Files

Status:

Version 1