Assessing Illumina technology for the high-throughput sequencing of bacteriophage genomes

Introduction

Viruses are the most abundant biological entities on the planet, having a ubiquitous distribution in all known biological niches. They play a crucial role in every environment they occupy, influencing the growth and metabolism of the hosts they infect (Ankrah et al., 2014; Mann et al., 2003; Meaden, Paszkiewicz & Koskella, 2015; Thompson et al., 2011), changing the chemical balance of the surrounding environment (Weitz & Wilhelm, 2012), affecting the diversity of the incumbent microbial community (Weitz et al., 2015), and driving microbial evolution through the transfer of genetic material (Lindell et al., 2005; Millard et al., 2004). In the marine environment, viruses are found at up to ~1 × 10⁸ ml⁻¹ of seawater, with an estimated ~4 × 10³⁰ viruses in the oceans (Suttle, 2005). Phages play an important role in biogeochemical cycling, diverting the flow of carbon to dissolved organic matter and particulate organic matter through the lysis of their bacterial hosts, thus influencing the amount of carbon that is sequestered to the deep ocean by the biological pump (Suttle, 2007). In soils, viral abundance is estimated to be 10⁸–10⁹ g⁻¹ (Williamson et al., 2013). The human gut is estimated to contain 10¹⁵ phages, although the exact role of phages in shaping the incumbent bacterial community structure and human health is unknown (Dalmasso, Hill & Ross, 2014). Unlike many other environments, such as the open ocean (Hurwitz & Sullivan, 2013), the majority of phages that have been found in gut metagenomes are temperate rather than obligately lytic (Howe et al., 2015). These phages possibly play an active role in human gut immunity and metabolism by forming part of the mucus-associated microbiome, where they may serve as a first line of defence against invading bacteria (Barr et al., 2013; Barr et al., 2015). In the mouse gut, phages have been implicated in providing a mechanism for multidrug resistance via reservoir genes (Modi et al., 2013). However, phages also hold promise as a possible alternative to antibiotics which are becoming ineffective with the increasing rise of antimicrobial resistance (Nobrega et al., 2015). In recent years, research of bacterial host immunity against phage infection has yielded one of the most promising genome-editing molecular techniques—the CRISPR-Cas system (Barrangou et al., 2007). Despite the importance of phages, there are relatively few isolates for which complete genome sequences are available. It is now over 30 years since the sequencing of the ssDNA phage ϕX174 in 1977 (Sanger et al., 1977), followed by the dsDNA phages lambda and T7 in the early 1980s (Dunn, Studier & Gottesman, 1983; Sanger et al., 1982).

The sequencing of these smaller phage genomes was completed many years before the first bacterial genome was published in 1995 (Fleischmann et al., 1995). In the last ten years, the number of bacterial and phage genomes has increased dramatically, coinciding with the decreasing cost-per-base of sequencing. Despite the smaller genome sizes of phages compared to their host, which should make them easier to sequence and assemble de novo, publicly available databases contain more finished bacterial than phage genomes. Whilst the numbers of finished bacterial genomes is slowly increasing (2,789 complete genomes within the European Nucleotide Archive—ENA, http://www.ebi.ac.uk/genomes/—at the time of writing), this is dwarfed by the tens of thousands of bacterial genomes that have been submitted as whole genome shotgun surveys (WGS) (e.g., ~40,000 Salmonella genomes alone). In comparison, there are 1,922 completely assembled phage genomes within the ENA, with no phage WGS assemblies in this database.

The vast majority of genes in viral metagenomes are not found in any of the currently sequenced phage genomes (Breitbart & Rohwer, 2005; Edwards & Rohwer, 2005; Hurwitz & Sullivan, 2013; Willner, Thurber & Rohwer, 2009). Increasing the number and diversity of sequenced phage isolates will improve our ability to assign unknown genes to a particular viral family. Combined with the additional information concerning the identity of the host they infect, the time of the year, and environmental conditions predominant at the location from which they were isolated, increasing the number of sequenced genomes will assist us in understanding the biology of a larger variety of novel phages. The recent sequencing of a small number of phage genomes has had a significant contribution towards understanding the bacterial metabolic processes in which they participate, highlighting the importance of further discovering the viral unknowns via sequencing of phage isolates (Chan et al., 2014; Chan et al., 2015; Kang et al., 2013; Roux et al., 2014; Sabehi et al., 2012; Zhao et al., 2013). The dearth of reference phage genomes infecting cells of a particular genera is highlighted by the estimation of 5,476 viral populations in the upper oceans, of which only 39 could be affiliated to cultured viruses (Brum et al., 2015). Re-sequencing of cultured phage genomes also provides insight into their genome evolution (Puxty et al., 2015). Altogether, it is clear that there exists a developing need for a high-throughput approach to the sequencing and re-sequencing of cultured phage genomes.

The limitations on the sequencing of phage isolates have previously been discussed in detail (Klumpp, Fouts & Sozhamannan, 2012). The largest single bottleneck will always be culturing of their host bacteria, which has been successfully circumvented by the use of metagenomics (Alavandi & Poornima, 2012; Blomstrom, 2011; Delwart, 2007; Edwards & Rohwer, 2005; Ge et al., 2013; Hurwitz & Sullivan, 2013; Kim, Whon & Bae, 2013). This development has greatly expanded our knowledge of phages, but there is ambiguity in identifying the hosts of phages found in metagenomes. Previously sequencing of phage genomes had been problematic due to cloning phage DNA for Sanger sequencing-based approaches, contaminating host DNA, and potential differences in %GC content relative to that of their hosts (Klumpp, Fouts & Sozhamannan, 2012), but these issues have been largely overcome by technological advances (Breitbart et al., 2003; Rohwer et al., 2001). Previous work, using 454 pyrosequencing and a column-based clean-up for DNA extraction has optimised methods for sequencing phages using this technology (Henn et al., 2010; Marine et al., 2011). Since then, 454 pyrosequencing technology has become largely obsolete with rapid advances in high-throughput sequencing technologies (see Loman et al. (2012), for a review). The Illumina MiSeq platform offers the potential to rapidly sequence hundreds of phage genomes. However, previous reports have suggested that Illumina data is of limited value for the de novo assembly of phage genomes (Klumpp, Fouts & Sozhamannan, 2012). Despite this, there are an increasing number of reports using Illumina technology for sequencing phage genomes (Carson et al., 2015; Malki et al., 2015; Smith et al., 2015; Tevdoradze et al., 2015). The resultant genomes have been assembled using a number of different assembly programs including SPAdes (Smith et al., 2015), Velvet (Cowley et al., 2015), and CLC Workbench (Carson et al., 2015) at a range of sequencing depths up to 18,000× coverage. Unlike their bacterial hosts, phages do not generally contain repetitive sequences such as gene duplications (e.g., rRNA operons), variable number tandem repeats or transposable elements that can prevent genome assembly. If these repetitive elements are longer than the library insert size a reliable assembly cannot be obtained (Magoc et al., 2013; Treangen & Salzberg, 2012). Thus, complete genome assembly should be possible using short read sequencing technologies, yet it is currently unknown what percentage of total phage isolates can be fully assembled, what the minimum coverage required for assembly is, and how likely they are to contain platform-specific assembly errors. Previous research has comprehensively evaluated the effect of both assembly program and depth of sequencing coverage for bacterial genomes (Magoc et al., 2013). However, the likelihood of a successful phage assembly, as well as the factors that may affect assembly, remain unknown. We aim to determine how the choice of assembler, depth of sequencing, and how multiplexing phages will influence the likelihood of a successful genome assembly.

Materials and Methods

Phage T4 DNA was purchased from Fluka. Using a CsCl purified stock of phage HP1, DNA was extracted. Cyanophage S-PM2d and S-RSM4 were cultured as previously described (Clokie et al., 2003). The remaining phage isolates were cultured on their respective hosts in King’s B medium (King, Ward & Raney, 1954) at 28 °C, purified with CHCl₃ and stored at 4 °C.

Bioinformatics Analysis

Phage genomes were downloaded from EBI in February 2013 and filtered to remove any genomes that contained unknown or ambiguous bases, resulting in 1826 genomes that were used for creating in silico datasets (see Table S1 for accession numbers). Simulated datasets of 2 × 300 bp paired-end reads were produced using ART Illumina 2.1.8 (Huang et al., 2012). Read quality profiles for input into ART were generated from a previous MiSeq run and produced error model profiles that are included as supplementary data. Insertion and deletion rates for in silico reads were based on the same read quality profiles. For initial comparison of assembly programs, 300 bp read sets were produced at 100× coverage of each genome, with a mean insert size of 300 bp. For further analysis, read sets were produced for a coverage of 20, 30, 40, 50 and 100×. Insert sizes of 300, 500 and 650 bp were used with selected datasets. Each set of reads was assembled with SPAdes v3.1 (Bankevich et al., 2012), Velvet v1.2.10 (Zerbino & Birney, 2008), and Ray v2.3.1 (Boisvert, Laviolette & Corbeil, 2010). For SPAdes the ‘–only-assembler’ parameter was used. For Velvet, a range of kmer values were tested with VelvetOptimiser using the parameters ‘-s 51 –e 199 –× 20 –cov_cutoff 4’. For Ray, the parameter ‘-k 61’ was used. All other parameters in each assembler were left at default settings. Genome assembly was assessed using QUAST (Gurevich et al., 2013) against the reference genome used for production of the reads. Genomes were defined as complete if a single contig assembled without error and covered > 97% of the reference sequence, as determined by QUAST (Gurevich et al., 2013). Partial assemblies did not meet the required 97% threshold and had no assembly errors. Misassemblies were defined as genomes that contained errors, as identified by QUAST analysis. To simulate the pooled libraries, genomes were randomly sampled with replacement from specific sets of phage genomes, either using all phage genomes or only those confined to a particular host (e.g., Pseudomonas, Mycobacterium, Synechococcus and Bacillus) prior to de novo genome assembly.

Results

To ascertain if the use of Illumina sequencing technologies is a suitable method for the high-throughput assembly of phage genomes, an artificial dataset was constructed from all high quality genomes (1826) in the EBI database at the time of analysis (http://www.ebi.ac.uk/genomes/phage.html). Genomes were assembled using three commonly used assembly programs (Velvet, SPAdes, and Ray) at a sequencing coverage of 100× to determine if the algorithm affects the successful assembly of phage genomes compared to the original reference genome. QUAST was used to assess the assembly parameters compared to the original assembly. Genomes were assessed on their completeness (% of genome on a single contig) and the number of assembly errors (Table S1). This resulted in 1800, 1545 and 1638 complete error free genomes, using SPAdes, Velvet, and Ray respectively (Table 1). A greater proportion of genomes were assembled with SPAdes, compared to both Velvet and Ray (Table 2). Velvet produced an increased proportion of assemblies that were misassembled, compared to the reference.

SPAdes provided the highest successful genome assembly at 100× coverage with an insert size of 300 bp; therefore, further parameters were tested using this assembler. Different levels of coverage (20, 30, 40, 50 and 100×) were then tested for each phage.

The ability to assemble genomes without misassembly is essential, if genome structure and synteny between phages is to be studied. The use of in silico reads from genomes that have previously been assembled allows the likelihood of misassemblies to be assessed. No misassemblies were found at 20–100× coverage, confirming that it is possible to correctly assemble phage genomes using SPAdes.

At 20× coverage, 98.0% of all phage genomes were assembled into a single contig (Fig. 1). At 30× coverage this increased to 98.5% (1800); increasing the coverage above 30×, provided no benefit to the assembly outcome (Fig. 1). Twenty-six phage genomes did not assemble completely even at 100× coverage, these phages had a range of genome sizes, %GC content, and different hosts (Table S2). Phages that did not assemble, had an average genome size of 121 kb, larger than phages that did assemble, with an average genome size of 66 kb. In order to try and successfully assemble these phage genomes into single contigs, we further investigated insert size and fold coverage in a range of combinations (Table S2). Insert size was increased to 650 bp and fold coverage to 1,000× (Table S2), resulting in the completion of one additional genome. The reason for the incomplete assembly in most instances is the failure to assemble the ends of the genome onto a single contig. For most genomes > 90% of the genome was assembled as a single contig, but this did not reach the 97% threshold (Table S2). Analysis of the genomes with ABACAS v1.03 (Assefa et al., 2009) against the known reference genome did not reveal any repeat regions that would prevent assembly.

These results indicate that it is possible to de novo assemble the genomes of the majority of phages correctly utilising current Illumina technology, under the condition that each phage sample was prepared as a separate sequence library. However, with a typical MiSeq run output being ~25 million paired-end reads, using a 96 multiplex strategy with dual indexes would result in an average coverage of > 1,000× of a 100 kb phage genome, or ~400× coverage using 384 indexes. A single MiSeq run has the capacity to sequence and assemble many times this number at 30× coverage. This capacity could be utilised through custom barcodes allowing the indexing of more samples. Such an approach carries with it the increased associated cost of library preparations using a larger number of indexes. An alternative approach is to pool phage genomes into a single library preparation. However, with de novo assembly, there is the possibility that the phages will be too similar to assemble into individual genomes post-sequencing.

In order to rule out the possibility of misassembly due to the similarity of different phages, two phage genomes were sampled (with replacement) from the 1820 phages previously assembled (Table S4), producing a model read set for each phage genome at an even coverage of 20×. This was repeated 960 times, thus simulating the output of 10 MiSeq runs. Combining two phages within a single library resulted in the complete assembly of both phage genomes in 98.8% of the samples (Table 3), when an even coverage was used with each phage genome. Misassembly occurred for eight genomes, when the phages that were combined were from the same or closely related bacterial hosts (Table 4).

To further test how co-sequencing genomes of similar phages will impact their assembly, additional simulations were carried out using combinations of sequences from phages isolated from the same bacterial host. Due to their numerical abundance in the current EBI dataset, phages infecting Mycobacterium (387 genomes) and Pseudomonas (142 genomes) were chosen for these simulations. Working on a basis of 96 libraries per MiSeq run, 96 pairs of Pseudomonas and Mycobacterium genomes were randomly sampled and combined, model reads were produced and genomes were assembled. Reads were combined at different sequencing ratios, with a minimum of 30× coverage at various ratios (1:1, 9:1, 1:9, and two randomly selected coverages). For both Mycobacterium and Pseudomonas phages, the effect of using different ratios of sequencing depth (compared to even coverage for both samples) was minimal; for Mycobacterium phages one extra genome was assembled at a 1:1 ratio compared to ratios of 9:1 or 1:9, while for Pseudomonas phages there was no difference in assembly outcome (Tables S5 and S6). For Mycobacterium phages it was possible to correctly assemble 166/192 genomes and partially assemble five genomes, with 21 genomes containing assembly errors (Table 5). For Pseudomonas phages, 171/192 phage genomes were correctly assembled, eight partially assembled and 13 contained assembly errors (Table 5). When compared to the random sampling of all phages, there was a clear increase in the proportion of phages that contained assembly errors when phages from the same host were combined. This finding presents a substantial issue when multiplexing novel phage genomes from the same host, since these assembly errors are more difficult to identify without prior knowledge of the genome sequence.

An alternative strategy for the high-throughput sequencing of multiple phages would be to isolate them from more than one host and combine them in a single library for sequencing. This was tested in silico, using Mycobacterium and Pseudomonas phages. This approach allowed the complete assembly of all 196 phage genomes without any assembly errors (Table 4). In order to test the limits of this multiplexing approach, this basic principle was extended in silico by production of read sets for mixtures of phages from three hosts (Mycobacterium, Pseudomonas and Synechococcus, see Table S7) and four hosts (Mycobacterium, Pseudomonas, Synechococcus and Bacillus, see Table S8). When combining three or four phages in a single sample, the proportion of phage genomes that were partially assembled increased with the complexity of the mixture, but crucially the number of misassemblies did not (Table 4).

In order to test the validity of the predictions reached via in silico experiments, we re-sequenced four known phages and compared the results of sequencing to the in silico predictions. We determined the ability to assemble each genome at a range of different sequencing depths by sub-sampling (Fig. 2). Phage T4 genome assembly required a threshold of 50× coverage and the median insert size was 350 bp (Fig. 3). In silico reads with an insert size of 300 bp allowed complete assembly of phage T4 at 20× coverage. Read mapping to the published T4 genome (accession AF158101) revealed 156 SNPs and 55 short INDELs (Tables S10 and S11). The large number of apparent mutations was initially surprising. However, the sequenced T4 DNA was commercially produced (Fluka) and the provenance of the T4 used in its production is unknown. The choice of reference will clearly affect the SNPs that are detected. There are currently seven complete phage T4 genomes (sensu lato); RB55 (KM607002), RB59 (KM607003), T4 (AF158101), T4 strain 147 (KJ477685), T4 strain GT7 (KJ477686), T4 strain wild (KJ477684) and T4T (HM137666), all of which have slightly different genome sizes. Each could serve as reference for SNP calling and produce differing results. For instance, T4 genes g10, g12 and g13 (Table S10) were all found to have at least one non-synonymous mutation, however, these would not have been called SNPs if T4T (HM137666) was used as reference. The identification of so many SNPs highlights the micro-diversity of what is referred to as phage sensu lato T4.

It was possible to assemble phage HP1 at 30× coverage, although coverage at the termini of the genome was strikingly lower, compared to the rest of the genome (Fig. 3). Additionally, a region with no coverage can be clearly identified at ~9,000 bp. Further investigation of the mapped reads identified a deleted region located at 9,287–9,406 bp when compared to the reference strain (accession: U24159). This causes an in-frame deletion of 39 amino acids from the sequence of the HP1p17 protein; whether the resulting protein is still functional remains unknown. Forty-three SNPs and six INDELs were identified in addition to the large deletion (Table S13). Twenty-nine SNPs caused non-synonymous mutations, eight were synonymous mutations and six were in intergenic regions (Table S12). Of the six INDELs that were identified, three occur in HP1p22, encoding a capsid scaffolding protein (Esposito et al., 1996) (Table S12). Although scaffolding proteins do not appear in the final capsid, they play an important role in its correct assembly. HP1p22 is a homologue of the GpO protein of P2, which is essential for the production of proheads (Lengyel et al., 1973). Therefore, despite the numerous INDELs and SNPs, this essential structural protein appears to still be functional, since this genome produced functional phages.

Cyanophage S-PM2d was assembled from two independent libraries from different DNA extractions, with a median insert size of 362 bp and 139 bp, producing assemblies designated S-PM2di and S-PM2dii respectively (Table 4). The larger insert library allowed assembly into a single contig at 50× coverage (Fig. 3), which was higher than the 20× coverage predicted by our in silico experiments. The 139 bp insert library required 100× fold coverage for complete genome assembly.

Whilst there was similar coverage at the termini compared to the rest of the genome, there was approximately double the number of reads that were mapped to the region 10,500–12,000 bp (Fig. 3). A similar pattern was observed for cyanophage S-RSM4, where a higher proportion of reads mapped to the region 1,37,800–1,39,000 bp. Assembly of cyanophage S-RSM4 was possible at 100× coverage (Fig. 2). Despite having similar %GC content and genome sizes (Table 6), the insert size of libraries for S-PM2i, S-PM2ii and S-RSM4 were very different, with median insert sizes of 362, 139 and 155 bp respectively (Table 6). Both S-PM2d and S-RMS4 contained no SNPs or small INDELs when compared to the published reference strains (accession LN828717 and FM207411).

The re-sequencing of four phages clearly demonstrates that complete de novo assembly of phage genomes from Illumina sequencing data is possible at coverage levels readily achievable on the MiSeq. For high-throughput phage genome sequencing, it will be necessary to remove other bottlenecks in the process. Standard NexteraXT library preparation only requires one ng of DNA for library preparation, as the protocol utilises transposons carrying adapter oligos to simultaneously fragment and ligate adapters in a single reaction (Marine et al., 2011). Researchers have optimised phage purification methods that use a column-based approach rather than time-consuming CsCl gradients (Henn et al., 2010). In this study, an alternative approach eliminated expensive columns and required only minimal phage lysate. This was achieved using a basic phenol:chloroform extraction (see methods) on a crude phage lysate.

At least one phage genome was assembled from each lysate (Table 6), with two genomes arising from the lysates of S-PM2dii, HC15b and VCM1. Assemblies used a minimum of ~1,40,000 reads, well in excess of the 30× coverage predicted using the in silico reads. After assembly, the total read pool was sub-sampled at random to produce datasets with coverage ranging from 20–100×. 30× coverage was sufficient to assemble complete genomes for seven of the nine novel phages, and the remaining two genomes assembled at 40 and 50× respectively (Fig. S1). In addition, 16 phage datasets produced using Illumina technology were extracted from the ENA (www.ebi.ac.uk/ena) and assembled at 20, 30, 40, 50 and 100× coverage (Fig. S2). All 16 of these genomes assembled as a single contig at 30× coverage.

Contaminating non-phage reads were identified by their failure to map to assembled phage genomes. The proportion of reads that could not be mapped to the phage genome varied from 0.01–47.69% (Table 6). In the case of S-PM2ii, somewhat surprisingly, two phage genomes were assembled from the lysate. The first was S-PM2d (~187 kb) and the second a genome of ~96.84 kb. The latter phage is of unknown origin, and not a cyanophage based on its lack of genetic similarity to known cyanophages. Furthermore, phylogenetic analysis of phoH, a common phage marker gene (Goldsmith et al., 2011), indicated that the closest phoH-containing relative is a phage infecting Caulobacter, not a cyanobacterium (Fig. S3). This phage is probably a result of a low-level contaminant in the stock of S-PM2d that was only detected due to the very high depth of initial sequencing coverage (> 2,000×). This contaminant has not been detected previously in S-PM2d which was recently purified by three rounds of plaque assay and sequenced to a high depth of coverage (Puxty et al., 2015).

Discussion

In silico predictions suggest that it is possible to correctly assemble the vast majority of phage genomes at 100× when using three common assembly programs. Using default setting with SPAdes, a greater number of genomes was assembled compared to using either Velvet or Ray. It is possible that further optimisation of parameters for both Ray and Velvet would increase the number of complete genomes. Further testing of coverage and insert size was carried out using SPAdes, but the results are likely to prove useful as a starting point for other assembly programs. Utilising SPAdes we found that increasing coverage above 30× resulted in a small increase of 0.5% in assembly success. However, 26 genomes were found not to assemble, even when a coverage of 1,000× or a large insert of 650 bp was used. Unlike their bacterial hosts, phages do not contain large repeats and multiple gene duplications that make bacteria recalcitrant to complete assembly (Ricker, Qian & Fulthorpe, 2012). For those genomes that were not assembled, it was the ends of the genome (as defined by the reference genome) that remained incomplete.

Whilst the use of NexteraXT has many benefits, the ability to consistently produce an insert size of > 500 bp for genomes of unknown size and %GC content is likely an unattainable target. As seen by the S-PM2d libraries tested in this study, the insert size was inconsistent, despite following the same protocol. Instead of trying to control the insert size, adjusting the genome coverage is a more attainable target, and one which can overcome very short insert sizes. Our experimental data showed that a higher coverage than predicted in silico was required to achieve complete genome assembly for some genomes. There are several possible explanations for this finding: (1) The reads that are produced in silico are completely random, unlike experimental data where some regions of the genome are clearly over-represented (Fig. 3); (2) The mean insert sizes of these experimental libraries was smaller than our in silico libraries (for example, in S-PM2, when the coverage required for assembly decreased from 100–40×; (Fig. 2); (3) A high proportion of non-phage DNA results in lower fold coverage than predicted, as was the case for S-PM2d and S-RSM4.

Using experimental read data from 13 phages in this study and 16 genomes retrieved from the ENA, we found that 82.75% (24/29) of these genomes could be assembled at 30× coverage, compared to the in silico data of 98.5% of genomes. Increasing phage sequencing throughput, by combining multiple phages within a single library preparation, would allow doubling of sequencing capacity at a marginal increase in cost for DNA extraction. The results of in silico predictions suggest that this would be a useful strategy if the phages were isolated from multiple hosts. This approach has previously been successfully used to assess the diversity of phage from the North Sea (Wichels et al., 1998). The recent development of a method of single-plaque phage sequencing (Kot et al., 2014), optimisations of the NexteraXT protocol to reduce the cost per library (Baym et al., 2015), and the results of this work all indicate it should be possible to take a high throughput approach to phage genome sequencing.

S-PM2d resequencing has confirmed, by chance, that it is possible to completely assemble the genomes of two phages from different hosts. Detection of the novel phage AM-2015 while sequencing S-PM2, was possible due to the extremely high sequencing depth used in this study, as it accounted for only ~13% of reads in the initial library. Similarly, in the case of phages VCM1 and HC15b, sequencing at high depth has revealed the presence of an additional phage in each of the lysates. All of these phages had undergone plaque purification prior to production of lysates. The presence of multiple phages within plaque purified isolates has previously been reported (Henn et al., 2010) and are known to be from the spontaneous release of a host prophage (Cowley et al., 2015). Without having complete host genomes for strains used in this study, it is not possible to determine if the presence of multiple phages is due to prophage release from the host.

Whilst it was possible to assemble multiple phage genomes from a single library in silico, simultaneous sequencing of multiple phages that are closely related could result in misassemblies that are extremely difficult to detect. Therefore, complete phage genomes assembled from metagenomes should be interpreted with care. We used SPAdes for genome assembly in this study, but acknowledge that it is not recommended for metagenome assembly. However, it seems likely that the misassembly of closely related phage genomes would result independent of the choice of assembler. When combining multiple phages from different hosts there was no increase in the number of misassembles, with a small increase in the percentage of incomplete genomes. Therefore, the multiplexing of phages from different hosts in a single library represents an efficient and cost effective way to increase the throughput of phage genome sequencing.

Conclusions

This work has demonstrated it is possible to assemble the vast majority of phage genomes using short read sequencing technologies in silico. Using 30× coverage was sufficient to assemble the majority of phage genomes both in silico and from experimental data. However, due to the uneven coverage produced by either the NexteraXT preparation or biological properties of the phage genome, a coverage of 100× would be recommended as a starting point to maximise the likelihood of successful assembly in a high-throughput manner. In silico analysis predicts that pooling two, three, or even four phages from different hosts into a single NexteraXT library preparation can increase the throughput of phage genome sequencing, without increasing misassemblies or additional costs for increased barcodes and library preparations. The current bottleneck in phage genomics is now clearly the ability to culture phages and isolate their DNA, rather than the capacity to sequence or assemble them.

Supplemental Information