Recombination within RNA viral genomes is a powerful driving force behind the evolution of RNA viruses and has been widely documented in bacterial, plant and animal viruses (for a review see 1 ).There are two possiblemodels for RNA recombination: non-replicative or breakage-rejoining;and copy-choice recombination. Breakage rejoining mechanisms are similar to that occurring during mRNA splicing and involve a catalytic ligation of two RNA molecules. During copy-choice recombination, RNA replication is halted mid-flow and the template RNA dissociates. The nascent RNA remains associated with its polymerase,however, andis able to re-prime replication either at a new position on the original template or on a new template strand. Copy-choice recombination is thought to be the most common mechanism for RNA recombination in many RNA viruses and is thought to primarily occur during the synthesis of the negative-sense RNA strand from the positivesense RNA genome.
RNA recombination may be advantageous to viruses through a variety of proposed mechanisms including the generation of viral escape mutants, the removal of deleterious mutations, the mating of viral genes and genomes or through the capturing of other viral or host genes.As a consequence, RNA recombination contributes to the extraordinary and rapid adaptability of RNA viruses and has been attributed to be the source of a number of recent outbreaks of, for example, dengue virus and echovirusas well as the emergence of entirely new viruses,including severe acute respiratory syndrome associated coronavirus (SARS-CoV). 11 Conversely, viral RNA recombination may also facilitate the stepwisegeneration of Defective-RNAs. 12 Defective-RNAs have lost their ability to independently encode functional viral proteins,but maintain the sequence motifs required for RNA replication by the viral RNA-dependent RNA polymerase (RdRp) and so are able to proliferate parasitically. Additionally, Defective-RNAs are often packaged into viral particles and so may contain the necessary motifs required for encapsidation. If a Defective-RNA is able to attenuate the replication of the full-length RNA genome through competition for access to the RdRp, then they are known as Defective-Interfering (DI) RNAs. 13 DI-RNAs arise during both persistent and acute infections most commonly in cell culturebut also occur during wild infections.
Flock House Virus (FHV) is anicosahedral T=3 non-enveloped virus containing a positivesense, single-stranded RNA genome. It is the prototypic member of the Nodaviridae family and provides an important model for the study of many pathogenic human viruses includingpoliovirus, human rhinovirus and hepatitis C virus (reviewed in). FHV naturally infects insects including Drosophila melanogaster as well as medically important insects such as mosquitoes (including Anopheles gambiae) and thetsetse fly. 20 The FHV genome consists oftwo positive-sense, single-stranded RNAs. RNA 1 (3.1 kb) encodes the RNAdependent RNA polymerase (RdRp) which can autonomously replicate the FHV genome. RNA 2 (1.4 kb) encodes the capsid precursor protein, α, 180 copies of which form the viral capsid.FHV packages one copy each of RNA 1 and RNA 2, but on average, nearly 2% of the packaged RNA is derived from transcripts of the host genome. 21 Additionally, virus-like particles (VLPs) of FHV can be made in cell culture by expressing the single capsid protein from a baculoviral vector in Sf21 cells. 22 VLPs spontaneously assemble into particles that are closely similarto native FHV, except that they package ribosomal RNAs, transcripts derived from the baculoviral expression vector and other cellular RNAs. 21 Next-Generation sequencing (NGS) has been used in multiple investigations to assess the polymorphism and mutation frequency within a viral genome, for example, in HIV, 23 human rhinovirus, 24 and foot-and-mouth disease virus. 25 With NGS technology it is possible to detect fusion events that may occur during RNA/DNA recombination and has been used to characterize mRNA splicing events 26 and to identify oncogenic chromosomal rearrangements.Here, we use NGS to investigate the sequence variation as a result of RNA recombination within the genome encapsidated by FHV. We profile with single nucleotide resolution hundreds of thousands of recombination events that have occurred during FHV genome replication and which have subsequently been encapsidated into virus particles.
FHV virions and VLPs were generated in cell culture and purifiedwith multiple sucrose gradients as is well established for structural and biochemical studies. 30 We thus obtained highly purifiedviral RNAfor RNAseq analysis by virtue of its encapsidation inside FHV virions and VLPs. We did not employ any additional RNA-selection steps, such as PCR or poly-A capture, which would yield a non-random population of RNA molecules. Moreover, the RNA encapsidated in FHV virions is directly relevant to virus infection as it is what may be delivered to another host cell during an infection.
The purified RNAswere fragmented and prepared for sequencing analysis using Illumina protocols for obtaining millions of single short reads.We thus obtained two raw datasets corresponding to the RNA encapsidated by FHV virions and VLPs -hereafter termed the FHV-RNAseq and VLP-RNAseq datasets. We applied a stringent quality filter to ourdatasets, as described in the materials and methods section, obtaining a total of 28.9M and 19.6M readsall exactly 95 nucleotides in lengthfor the FHV-RNAseq and VLP-RNAseq datasets respectively.
To establish the identity of the RNA reads in the datasets, we mapped single reads using the Bowtie aligner 31 to genes already known to be present in either the FHV virions or VLPs.We have previously thoroughly characterized the contents of FHV virions and VLPs, and so used these sequences as references here. 21 Reads were mapped in an end-to-end and un-gapped manner. The number of mapped reads forboth datasetsis shown in Table 1 . FHV virions packaged roughly equal stoichiometric amounts of FHV RNA 1 and RNA 2.Out of the 24.6 M reads that aligned to the FHV genome, only 4 of these mapped to the negativesense RNA transcript.In addition to the FHV genome, we also detected some encapsidated host RNA, as we have recently reported. 21 VLPs primarily packagedhost RNA, most of which consisted of ribosomal RNA.VLPs also packaged transcripts from the baculoviral expression vector, which included the FHV RNA 2 gene encoding the FHV capsid protein.
We also detected a small amount of FHV RNA 1 in our VLP-RNAseq (1,103 reads), indicating that there was a low-level persistent infection with FHV in our Sf21 cells, as has also been reported for a number other cell lines.Consequently, we would also expect approximately 300 of the FHV RNA 2 reads from the VLP-RNAseq dataset to be present as a result of FHV replication (judging from the stoichiometric ratios of FHV RNA 1 and RNA 2 seen in the FHV-RNAseq dataset), which is well below the total number observed in the VLP-RNAseq dataset (1.46M). Note that the FHV RdRp would not be able to replicate the FHV RNA 2 derived from the baculovirus vector as this is cloned without the appropriate 5' and 3' UTRs required for viral replication. 6 To illustrate the detection of recombination, we mapped the reads to two specific Defective-RNAs(R2D675 and R1D1626) that are known to arise during passaging of FHV in cell culture.These Defective-RNAs are formed by recombination events that occur between different portions of the FHV genome as illustrated in Figure 1 . Many of these result in large deletions of the RNA genome while others result in insertions or duplications (e.g. in R2D675, nucleotides 102-154 are duplicated and re-inserted near the 3' terminus). When these recombination events occur, a junction is generated with a novel nucleic acid sequence. From the FHV-RNAseq dataset, 15 single reads mapped to R1D1626, and 101 mapped to R2D675, all of whichmapped to junction sites that characterize the respective Defective-RNA. The locations of the mapped reads are illustrated in Figure 1 . As expected, no reads aligned to the portions of the Defective-RNAs that maintain a wild-type sequence as these reads were already aligned to the fully intact FHV Genome. However, not all junction events were detected, indicating that the precise character of the Defective-RNAs generated in our sample was different to that previously identified. The mapping of reads over junctions present in known Defective-RNAs is important to illustrate how individual RNA recombination events within the FHV genome can be detected.
We extended our analysis to detect novel recombination events occurring within the FHV genome.To achieve this, we generated a referencepseudo-library containing millions of short reference sequences (≤150 nucleotides in length) each corresponding topotentialrecombination events within the FHV genome.Thepseudo-library was generated using a Python script whereby every 75 base sequence from the FHV genome was appended to every other 75 base sequence. We also allowed shorter sequences to be generated for reference sequences corresponding torecombination events occurring near the edges (≥20 bp) of the FHV genes. The reference pseudo-library contained 19,965,801 sequences describingjunctions within FHV RNA 1 and FHV RNA 2 as well as between these two. For each reference, the first 75 nucleotides on the 5' side of the junction site is defined as the 5'strand, while the 3' nucleotides form the 3'strand.Thus, by mapping 95-nucleotide reads to reference sequences 150 nucleotides in length where the recombination junction is placed between bases 75 and 76,a minimum of 20 nucleotides from the single readmust align to boththe 5' and 3'strands.Junctions appearing in either 20 bp extremity of the single reads would not bedetected. Consequently, there are 56 of a possible 94 'cutting' sitesin the 95 bp single reads where junctions maybe detected.While this will result in an under-reporting of all possible recombination events, this prevents the possibility of mapping reads with too few nucleotides on one side of a junction to unambiguously assign its identity.
The single reads from both the FHV-RNAseq and VLP-RNAseq datasets that did not align to the wild-type full length reference sequences were next mapped to the junction reference library ( Table 1) . Every read that mapped to the junction reference library mapped to the positive-sense strand.In the FHV-RNAseq dataset, we detected a wide variety ofjunctions within each of the FHV genes. The majority of these were foundwithin RNA 1, accounting for approximately 4% of the single reads from the FHV-RNAseq dataset and primarily correspondedto large deletions. A smaller number(~0.1%)were found in RNA 2, and again, primarily corresponded to deletion events.Junctionswere observed that result in effective insertions or duplications of the FHV genes(where the 5'strand is down-stream of the 3'strand).Junctions were also detected betweenFHVRNA 1 and RNA 2 indicating that RNA recombination between non-homologous templateshas occurred, although these events are less frequent than those within RNA 1 or RNA 2 despite there being a greater number of possiblerecombination events that can be mapped.
These data are represented as heat maps in Figure 2 . Here, the y-axis describes the last nucleotide of the 5' strand and the x-axis describes the first nucleotide of the 3' strand. The number of reads that map to each particular junction is indicated with a color bar. The heat maps illustrate the widearray of junctionsdetected, reflecting the diversity of the packaged RNA as a result of recombination.
Prominent horizontal and vertical striations are visible in the heat maps in Figure 2 . Interestingly, the locations of the striations are maintained regardless of the identity of the recombination partner; for example, the horizontal striation occurring between nt 710-730 in RNA 2 is present regardless of whether RNA 2 or RNA 1 provides the 3' strands. These striations may therefore indicate a sequence-dependent preference in the selection of recombination sites. It is interesting to note that some of these striations pass through locations of previously characterized Defective-RNAs as well as locations of high-frequency events detected here. The positions of the horizontal and vertical striations do not correlate (Pearson Coefficients are -0.0035 and -0.0019 for RNA 1 and RNA 2 recombinations respectively). This indicates that the sequence-dependent selection of 5' strands differs from that of 3' strands. However, no apparent nucleotide preference at either strand could be detected, other than a weak under-representation of guanines at the junction sites ( Figure  S1 ). Additionally, there was no average sequence identity detected between nucleotides upstream of the junction site in the 5' strand and nucleotides upstream of 3' strand, as has been reported for Defective-RNAs, for example, in dengue virus 16 . The selection of recombination sites may therefore have a more complex origin, for example, in the character of the local RNA secondary structure. 6 The positions in the heat map that correspond to mapping to the Defective-RNAillustrated in Figure 1 are shown in Figure 2B and C. The precise sites of the junctions corresponding to the previously identified Defective-RNAs are marked with blue cross-hairsand are relatively poorly represented (the reads indicated under the cross-hairs in the heat map are the same reads illustrated in Figure 1 ). However, we detected twohigh frequency junctions nearby that join nt 249-517, and nt 730-1,229 (8,544 and 1,080 mapped reads respectively -these can be seen in Figure 2B insetsf and g respectively) and are two highest frequency events detected within RNA 2. This indicates that a different yet similar population of Defective-RNAs was favored in our virus sample.
We detected many high frequency events in FHV RNA 1 (58 unique junctions with >1,000 mapped reads). Only one of these resulted in an effective insertion event (from nt 1,605 back to nt 1,087). The remainder resulted in deletions ranging from 13 to 1,063 nucleotides in length. Figure 3 shows two regions of RNA1 to RNA1 recombination that are enriched with high-frequency events (349,194 reads in Figure 3A and 224,347 reads in Figure 3B ). Interestingly, almost the entirety of these events resulted in deletions exactly 3n nucleotides in length (99% of the reads shown). As a result, the ORF was maintained in each of these cases suggesting that the RNA arising from these recombination eventshas been selected by virtue of their ability to produce a viable protein product. This may be either because RNA encapsidation is coupled to translation(as has been demonstrated for FHV 35 ) or because a functional, yet truncated, form of the viral RdRp is being selected.
The high-frequency of these recombination events suggests that they are present due to successive rounds of replication rather than independent instances of RNA-recombination. This is supported by the fact that the distribution of frequencies with which individual recombination events were detected is heavy-tailed, reminiscent of a power-law distribution (Figure 3C-F) . This indicates that there is an over representation of a small number of unique recombination events. Heavy-tailed distributions can arise when events are initially randomly generated, but some selected events are favorably duplicated, 36 as has been observed for a number of replicable components of eukaryotic and prokaryotic genomes. 37 Such a scenario is what we would expect during the generation of Defective RNAs, which are initially generated by stochastic RNA recombination but are subsequently highly (sometimes competitively) replicated by viral RdRps. 12
Many of the recombination events that we detected above occurred within close proximity to one another and it is likely thatmany of these were present on the same originalDefective RNA molecule.Consequently, single reads would be present in our dataset that will span two recombination events and thus would not be able to map to our pseudo-library containing reference sequence with only single recombination events. To address this, we designed a second pseudo-library of reference sequences that contained combinations of twopreviously detected recombination events that occurred within close proximity. As before, we designed this library to enforce the single reads to map with at least 20 nt on the 5' and 3' sides of the recombination events. Consequently, the maximum distance allowed between junction sites was 55 nt (95 nt from a single read minus two 20 nt 'seeds') and we allowed a minimum of 5 nt. This second pseudo-library contained 19'379'090 reference sequences ranging from 95 to 145 nt in length.From the 1,990,779 reads that remained unaligned ( Table 1) , we mapped an additional 80,374 reads to RNA 1 recombinations and 58 reads to RNA 2 recombinations. As each reference contained two recombination junctions, this corresponds to an extra 160,748 and 116 junction sites detected in RNA 1 and RNA 2 respectively.
We also evaluated the VLP-RNAseq dataset for the presence of recombination events. The VLPs were generated by expressing the capsid protein (FHV RNA 2) from a baculoviral expression vector.In the absence of viral RNA replication, VLPs package host RNA and RNA transcripts from the baculoviral expression vector ( Table 1) ,which are transcribedby the host DNA-dependent RNA polymerases. We would not therefore expect any junctions to be present in the VLP-RNAseq dataset as a result of recombination during viral RNA replication. The VLP-RNAseq dataset was generated using the same procedure as for the FHV-RNAseq dataset.Consequently, by searching for junctions within the VLP-RNAseq dataset we can estimate the amount ofrecombinationthat has occurredduringPCR steps in the cDNA library preparationused for RNAseq, as is a recognized artifact. 38 From the VLP-RNAseq dataset, 82 reads mapped to RNA1-RNA1 junctions(note that RNA 1 was present in very small quantities due to a potential low-level persistent infection of Sf21 cells with FHV) and 2,168 reads mapped to RNA2 to RNA2 junctions (Figure 4 and Table 1 ). A prominent diagonal striation is visiblespanning theheatmap of RNA 2 to RNA 2 junctions.These events correspond to very short insertions and deletions. A histogram of the lengths of insertions and deletions that occurs due toeach recombination event ( Figure 5) shows that insertions and deletions 5 nucleotides or shorter (known as MicroInDels) areabundant in the VLP-RNAseq dataset as well as in the FHV-RNAseq dataset. This indicates that their formation was not unique to the FHV viral polymerase, but was most likely to have arisen during the amplification of the cDNA library used for sequencing. DNA polymerases are known to accrue MicroInDels and have been reported to be abundant in other NGS datasets Consequently, we exclude all MicroInDels when counting and comparing recombination events between the FHV-RNAseq and VLP-RNAseq datasets, as indicated in Table 1 . Excluding these MicroInDels, only 356 reads from VLP-RNAseq dataset mapped to recombination events in FHV RNA 2 (none were detected in FHV RNA 1). In contrast, we detected 26,958 junctions in FHV RNA 2 within the FHV-RNAseq dataset.
In addition to the FHV genome, we also made reference libraries for junctions occurring within two other highly abundant genes that were present in the VLP-RNAseq dataset: ORF '1629' from the baculovirus expression vector and a portion of the 18S rRNA. Together with FHV RNA 2, these transcripts made up a similar proportion of the VLP-RNAseq dataset (4.19M of 19.6M reads = 21.3%) as FHV RNA 2 alone did in the FHV-RNAseq dataset (6.75M of 28.9M reads = 23.4%). They thus provide a suitable control against which to compare the FHV-RNAseq dataset. We detected 1,220 reads that mapped to recombination events occurring within these genes and between these genes (Table S1) thus giving a total of 1,576 artifactual recombination events per 4.19M mapped reads. This would be equivalent to 2,539 reads per 6.75M mapped reads.Excluding the double recombinations, we detected 26,958 recombination events within FHV RNA 2 (6.75M mapped reads) and 1,167,929 reads within FHV RNA 1 (17.9M mapped reads) and a total of 33,356 between these two, which is clearly in great excess of the artifactual recombination detected in the VLP-RNAseq dataset. We can therefore be confident that the background recombination noise in our FHV-RNAseq dataset is low and that the recombination events observed in the FHV-RNAseq dataset primarily reflect the activity of theviral RNAdependent RNA polymerase.
The coverage of reads across the wild type FHV genome was not constant (Figure S2) , as is common in RNAseq studies owing to PCR-mediated bias.Consequently, the number of detected junctions over FHV RNA 1 and 2 were normalized by dividing the number of 5' strand recombination events or 3' strand recombination events by the number of wild type reads that mapped at each nucleotide position to obtain the frequencies of recombination events. From this we can make an estimate of the average recombination frequency across the FHV genome ( Table 2) . Thiswill reflect both the frequency of individually generated recombination events as well as the replication of junctions that are found in the Defective-RNAs. However, the observed frequencies provided a metric with which to evaluate the amount of genetic variety in thevirus sample.As these frequencies are likely to be inflated at the extremities of each gene due to the low coverage of wild type reads in these regions (in particular at the 3' terminus) (Figure S2) , we also show the rates over just the ORFs of each gene.These values indicate that RNA recombination is an abundant source of mutagenesis and is comparable in magnitude to that of mismatch mutation ( Table 2) .
As we are using RNAseq with single short reads, we are unable to detect recombination events that occur between two homologous templates but at identical sites as this would result in the conservation of the local nucleic acid sequence. However, the high rate of RNA recombination that we do detect suggests that such 'silent recombination' will also be abundant. Indeed, such a process would be highly important by allowing for the mating of homologous templates, potentially removing deleterious mutations or by allowing the reshuffling of advantageous mutations that have occurred on separate RNA molecules.
Next-Generation sequencing (NGS) has proven itself to be a valuable tool in assessing the mutational landscape of a viral genome. NGS has been used to map the positions of singlenucleotide polymorphisms in viral populationand to measure the frequency of mismatch mutation.Both of these are source of considerable diversity within the 'genetic cloud' of a viral genome and are used to characterize the quasi-species present in a virus sample. Here, we have laid out an approach that extends these capabilities to include RNA recombination. By mapping the position and frequency of every possible junction within the genome of FHV, we present a highly-detailed and complex landscape of the numerous recombination events that occur during viral RNA replication.
The strategy laid out in this paper could equally be applied to a wide range of virus samples, including DNA viruses, and could add to our understanding of their diversity and evolution. The frequency of RNA recombination is known to vary widely between viral species, even among different positive-sense RNA virus, and could be assessed using NGS. The frequency of RNA recombination could also be compared between different preparations of the same virus, e.g. during the course of an infection, or when amplified in different host cells for viruses that have a broad host-range (e.g. dengue virus) or when viruses are exposed to antiviral therapies (e.g. ribavirin treatment for HCV). Additionally, it would also be possible to search for recombination between different viral species during co-infections. This would be important for understanding the role of RNA recombination in the evolution of new viruses.
Our approach could also be used to discover and characterize novel Defective and Defective-Interfering (DI) RNAs potentially present in a variety of infectious viruses.The generation of DI-RNAs has been proposed to be a critical stage in the transition of acute to chronic viralinfections 43 and DI-RNAs have been found in patients persistently infected with measles virus, 44 dengue virus 45 and hepatitis C virus. 46 Characterizing DI-RNAs present even in very low titers may improve our understanding ofviral infections and helpidentify variations of such elements between individuals or host organisms or during the progression of a viral infection.Additionally, characterizing what portions of the virus genome are present in the DI-RNAs will help us understand which components of the genome are necessary for replication by the viral polymerases.Similarly, by analyzing the nucleic acids packaged inside viruses, we may potentially find which components of the genome are required for successful encapsidation. Finally, therapeutic applications could be envisioned as discovering the identity of DI-RNAs may allowfor their exploitation for treatment or prevention of acute viral infections. Such applications have been demonstrated, for example, in the form of deliverable vaccines 47 or through the transgenic expression of DI-RNA-like molecules. 48
For authentic Flock House Virus production, cultured S2 cells were grown in Schneider's media (Sigma) containing 15 % fetal bovine serum (Gibco) using standard laboratory procedures. Cells were concentrated to 4 × 10 7 cells per ml, infected with Flock House Virus at an MOI of 1 and rocked for 1 hour at room temperature. The cells were then diluted with Schneider's insect media to a final concentration of 8 × 10 6 cells per ml and incubated at 27°C in a rotary shaker. Cells were harvested two days post-infection. For the production of Virus-Like-Particles of Flock House Virus, Sf21 cells were cultured in TC-100 media (Invitrogen) supplemented with 10 % fetal bovine serum (Gibco) using standard laboratory procedures. FHV RNA 2 was expressed from the pBacPAK9 baculovirus vector as described. 22 Cells were harvested three days post-transformation.
Authentic virions and VLPs were purified using a series of sucrose gradients as is well established 30 in 50 mM Hepes pH 7.0. Clarified cell lysates were spun at 40,000 RPM for 2.5 hours onto a 30 % sucrose cushion. The pellet was resuspended and then applied to a 10-40 % sucrose gradient and spun at 40,000 RPM for 1.5 hours. Fractions from the sucrose gradient were removed and analyzed by SDS-PAGE. Fractions containing only viral capsid proteins were pooled. This sample was then supplemented with 10 X DNase I reaction buffer (NEB) and 20 Units of DNase I (NEB) and 0.5 μg of RNase A (Roche) and incubated at room temperature for 2 hours to remove any non-encapsidated co-purified DNA or RNA. The samples were then extensively washed with 50 mM Hepes pH 7.0 on a 100 kDa MWCO centrifugal concentrator. This sample was then applied to the top of a second 10-40 % sucrose gradient and spun at 40,000 RPM for 1.5 hours. Again, fractions from the sucrose gradient were removed and analyzed by SDS-PAGE. Fractions containing only viral capsid proteins were pooled and extensively washed on a 100 kDa MWCO centrifugal concentrator. After this extensive purification, no cellular proteins could be detected by coomassie stain on an SDS-PAGE gel and no RNA or DNA could be detected on a native agarose gel when loading 3 μg of virus.
Purified FHV or VLP particles were disrupted at room temperature by incubation in 0.1 % SDS and 0.1 M NaCl for 15 minutes. RNA was extracted from the disrupted particles using an equal volume of acid phenol followed by three washes with 100 % chloroform. RNA was then ethanol precipitated in the presence of 100 mM Sodium Acetate pH 5.3. RNA pellets were washed in 70 % ethanol, dried, and resuspended in pure water.
0.4 μg of RNA was prepared for Next-Generation Sequencing using a modified version of the Illumina protocol (http://www.illumina.com/applications/sequencing/ rna.ilmn#strand_specific_rna_seq) where 12 cycles of PCR were performed and standard Truseq adapters and Truseq barcoded primers were used. A final size selection was performed by native agarose gel electrophoresis to yield a library of inserts approximately 200 bases in length suitable for 100 base single read sequencing. The library was extracted from the agarose gel using standard oligo purification columns. The prepared library was then loaded onto an Illumina HiSeq v3 single-read flowcell, standard cluster generation was performed on a Cbot and sequencing was performed for 100 bases of the insert and 7 bases of the index read using standard HiSeq sequencing reagents on an Illumina HiSeq 2000 instrument. Reads were processed using CASAVA 1.8.2 and demultiplexed based on index sequences.
Reads containing any fragment of the 3' Truseq adapter were detected and trimmed using cutadapt (http://code.google.com/p/cutadapt/) 49 with default settings. Reads containing any base with a PHRED score < 20 were discarded using the FASTXtoolkit (http:// hannonlab.cshl.edu/fastx_toolkit/). The quality of the reads in the dataset was assessed using the FastQC package (http://www.bioinformatics.bbsrc.ac.uk/projects/fastqc/). This revealed a poor average base-calling quality in the final 5 nucleotides of each read. Consequently, each read was trimmed down from the 3' end to a total of 95 nucleotides in length. Shorter reads were discarded.
The Drosophila melanogaster reference genome r5.22 was downloaded from the FlyBase repository (http://flybase.org/) and the mRNA refSeq library was obtained from the UCSC genome browser website (http://genome.ucsc.edu/). ESTs from Spodoptera frugiperda cell lines were downloaded from the database: "Spodobase" (http://bioweb.ensam.inra.fr/ spodobase/). 50 Sequences for the FHV RNA 2 (NC_004144), FHV RNA 1 (NC_004146), the Attacusricini 45s rDNA (AF463459), the Baculovirus genome (NC_001623) and
Defective-RNAs (GU393238 and GU393241) were obtained from NCBI. Reads were aligned to the host genome reference using the Bowtie alignment package version 0.12.7 (http://bowtie-bio.sourceforge.net/index.shtml) 31 in -v mode, tolerating up to 3 mismatched nucleotides per 95 base read. Alignment files were processed using SAMtools (http:// samtools.sourceforge.net/) 51 and alignments were visualized and inspected using Tablet. 52 Reads were mapped to the FHV genome using Bowtie parameters: -v 2 --best. Junctions were detected by alignment of the remaining reads using Bowtie parameters -v 2 --best to a library of sequences corresponding to all possible recombination events in the FHV genome as described in the main text. Reads that mapped with mismatches and that mapped to the edges of the reference sequences were removed from the alignment (from the .sam file). This is because the mismatch tolerance can allow a read to map with fewer than the required 20 nucleotides at either the donor or acceptor strand by claiming mismatching at the junction site of an adjacent but incorrect reference.
Refer to Web version on PubMed Central for supplementary material.
• Diversity in RNA packaged by Flock House Virus was assessed by deepsequencing
• RNA recombination was quantitatively profiled with single-nucleotide resolution
• RNA recombination is prolific, with a frequency comparable to mismatch mutation
• Previously known, as well as novel, Defective RNAs are identified • Our strategy may be applied to detect recombination in a wide range of viruses Defective-RNAs that have been previously characterized in FHV infections are formed by recombination withinA) FHV RNA 1 to form R1D1626 and B) FHV RNA 2 to form R2D675. These recombination events generate unique junction sequences,labeleda-i for reference in Figure 2 . Reads from the FHV-RNAseq dataset were mapped to these sequences and are indicated underneath the schematic for each Defective-RNA to illustrate their mapping over the recombination junction sites. This reveals clusters of high-frequency events that potentially represent Defectively RNAs.Log-log plots of the ranked frequencies of unique recombination events between C) FHV RNA1 to FHV RNA 1, D) FHV RNA2 to FHV RNA 2,E) FHV RNA2 to FHV RNA 1and F) FHV RNA1 to FHV RNA 2 indicate that their distribution is heavy-tailed. A heat map similar to those illustrated in Figure 2 of the junctions detected in the VLPRNAseq dataset demonstrate that the majority of events correspond to MicroInDels (N = 1,812), evident as the strong diagonal striation. Other artifactual junctions were also detected, but with low frequency (N = 356). 