Lab-infected mosquito studies (Bishop-Lilly et al., 2010; Hall-Mendelin et al., 2013) suggest the feasibility of using SMS to identify both viruses and other microorganisms (such as bacteria and fungi) in wild populations of mosquitoes, however, to our knowledge, this has never been done. Because symbiont mediated transmission disruption is an emerging tool in controlling vector-borne diseases (reviewed in Cirimotich et al., 2011; Weiss and Aksoy, 2011), SMS of wild mosquitoes may be particularly informative because it simultaneously identifies both viruses and other symbionts and could uncover any correlation between the two in the same host. This study represents a proof of concept in characterizing microbiota using a taxonomically broad approach that could ultimately prove useful in exposing significant novel interactions between microbes. In particular, we use SMS to characterize the microbiota of seven individual mosquitoes of three different species that were collected from a variety of natural habitats in northern California, USA, which is an area where *West Nile virus* has been reported. Using RNA-based SMS on the Illumina platform, we identify sequences related to viruses, bacteria, and fungi in each individual. Furthermore, we were able to verify mosquito species identities using SMS data alone. This work exemplifies the utility of SMS to study the natural microbiota of mosquitoes and we hope it prompts future research in this area.

### Materials and Methods

### Collection and Identification

Several hundred mosquitoes were collected from seven locations in northern California (Pepperwood Preserve, Bolinas, Stinson Beach, San Rafael, Mill Valley, San Francisco, and San Mateo) during March–November 2013 (**Figure 1**). Locations were chosen to represent a range of habitats, from sylvatic/wild (e.g., Pepperwood Preserve) to urban (e.g., San Francisco), (**Table 1**). Collections occurred under the permit and permission

FIGURE 1 | Collection locations and a summary of the microbial sequences identified in each sample. Locations A through G correspond to the samples as listed in Table 1. A summary of the viral, bacterial, and fungal sequences identified in each sample is

agreements of the Marin/Sonoma Mosquito and Vector Control District, or on private lands with the owners' permission.

We employed several methods to collect mosquitoes, including hand nets, gravid traps baited with hay-infused water, or ZumbaTM traps baited with carbon dioxide gas (CO2), heat included. Details for each microbial group can be found in Tables 2–4. The map was created using Google Earth Pro, with data from the U.S. Navy, LDEO-Columbia, NSF, NOAA, SIO, NGA, MBARI, and GEBCO.

packs, and SkinLureTM. All samples were frozen within 24–48 h of being trapped and stored at −80◦C. We identified mosquitoes morphologically using a dissecting microscope and key (Bohart and Washino, 1978), retaining voucher specimens for each species and trapping event. In several cases, morphological

#### TABLE 1 | Sample details.

#### TABLE 2 | Viral sequences identified in northern California mosquitoes.

a*In amino acids.*

b*To nearest match within the custom blastp database (see Section "Materials and Methods" and supplementary data at http://dx.doi.org/10.6084/m9.figshare.1247641).* c*As determined by nearest blastp match and phylogenetic reconstruction.*

d*As determined using by mapping the quality-checked reads back to the assembled contigs using Bowtie 2. Depth is the average coverage across the entire length of the contig. Coverage maps for each contig are available at http://dx.doi.org/10.6084/m9.figshare.1247641* e*RNA-dependent RNA polymerase (RdRp).*

f *Proline-alanine rich protein (PArp).*

identification was verified by sequence identity for cytochrome oxidase I (COI) using a leg from the voucher specimen.

### DNA Extraction, Library Preparation, and Sequencing

From the several hundred mosquitoes collected, one representative individual female from each of the seven different collection sites was selected for SMS sequencing (**Table 1**). Prior to processing, whole mosquitoes were washed in 70% ethanol, distilled water, and phosphatebuffered saline (PBS) solution to remove external microbes. While we did not explicitly test the final wash for complete removal of external microbes, a similar wash protocol for *Drosophila* was found to be sufficient for this purpose (Chandler et al., 2011). Washed samples were then individually homogenized in PBS with steel Lysing Matrix I beads (MP Biomedicals).

Samples were prepared for sequencing by first undergoing total RNA isolation, rRNA subtraction, reverse transcription, and random amplification into cDNA libraries. For RNA extraction, we used the MasterPureTM Complete DNA and RNA Purification kit according to the manufacturer's protocols. After DNase treatment, pellets were resuspended in 30 µl of TE buffer and 1 µl of RiboGuardTM was added to each tube. Next, we performed ribosomal subtraction using Ribo-ZeroTM Gold (magnetic beads, human/rat/mouse kit). While we did not design our experiment to explicitly test the efficacy of Ribo-ZeroTM Gold on the removal of mosquito ribosomal TABLE 3 | Bacterial sequences associated with northern California mosquitoes.

a*RDP Family level confidence.*

b*RDP Genus level confidence.*

c*Manually confirmed to be Wolbachia by querying to the NCBI non-redundant database.*

d*Multiple contigs from the same sample match to the same bacterial genus.*

e*Only the bacterial portion of this chimeric contig was used.*

RNA (rRNA), we note that the final percentage of mosquito 28S and 18S rRNA per library was, on average, 29.5% and 2.5%, respectively (data not shown). This was determined using the default settings of Bowtie 2 (Langmead et al., 2009) to map the quality checked reads (see below) to any contigs with a closest match to insect large subunit (LSU) and small subunit (SSU) rRNA.

Samples then underwent random reverse transcription followed by random amplification. Samples were tagged with a molecular ID tag (MID) and pooled equimolar for a single flow cell of sequencing. Sequencing was done on the Illumina HiSeq 2000 at HudsonAlpha Institute for Biotechnology, using an insert size of 100 bp paired-end reads. No water or blank sample was included as a negative control [although this is recommended for future sequencing runs (Salter et al., 2014)]. Raw, unprocessed sequencing reads are available through the NCBI Short Read Database as part of BioProject PRJNA269777.

### Sequence Processing and Taxonomy Assignment

An average of 35 million reads were produced per library (details for each library are available in **Table 1**). Sequences were quality checked and assembled using the A5 pipeline with the metagenome flag (Tritt et al., 2012). A5 combines sequence quality control, adapter trimming, and contig assembly. The SGA software package removes low quality reads and corrects sequencing errors1 and Tagdust removes sequencing adapter contamination (Lassmann et al., 2009). Cleaned sequences are used to build contigs with the IDBA-UD assembler (Peng et al., 2012).

To identify any viral sequences associated with these mosquitoes, contigs were translated into all six frames and any open reading frames (ORFs) longer than 100 amino acids in length were further examined. These ORFs were queried against a custom blast database using the blastp algorithm. We used

1https://github.com/jts/sga/tree/master/src

TABLE 4 | Fungal sequences associated with northern California mosquitoes.

a*RDP Family level confidence.*

b*RDP Genus level confidence.*

c*Manually confirmed to belong to the family Cordycipitaceae by querying to the NCBI non-redundant database.* d*Multiple contigs from the same sample match to the same fungal genus.*

an *e*-value cutoff of 1 × 10−3, which is approximately equivalent, for our database, to 30% similarity over 100 amino acids and 20% similarity over 500 amino acids. This custom database contains the entire NCBI non-redundant protein database, the PhAnToMe phage protein database (Aziz et al., 2012), the NCBI Viral RefSeq Database, and the *Aedes aegypti*, *Anopheles gambiae*, and *Culex quinquefasciatus* protein databases (Megy et al., 2012). The rationale for combining databases was to ensure that viral and phage proteins were present, while simultaneously reducing the false positive rate by the inclusion of all possible non-viral or non-phage proteins. While none of the mosquitoes examined in this study have their complete genomes sequenced, we note that this would not reduce our ability to detect viral sequences; rather this simply increases the false positive rate. All ORFs with a closest match to "virus" or "phage" were manually queried to the NCBI website to confirm their identity. Confirmed viral sequences are available on NCBI through the GenBank accession numbers KP642114 to KP642128.

To identify any bacterial and fungal sequences in the datasets, the non-translated contigs were queried using the blastn algorithm (Altschul et al., 1990) to the SILVA SSU and LSU Reference Databases (Release 111; Quast et al., 2013). Contigs with a closest match to "Bacteria" and "Fungi" when queried to the SSU and LSU databases, respectively, and were longer than the 300 bp, were then submitted to Ribosomal Database Project's (RDP) Classifier for taxonomic assignment (Wang et al., 2007). Contigs above a 90% genus (SSU) or 70% family (LSU) confidence cutoff were considered reliable hits (**Tables 3** and **4**). All intermediate files, include blast results and fasta files of significant hits are available at http://dx.doi.org/10.6084/m9.figshare.1247641.

To identify the mosquito host taxa, the non-translated contigs were queried to a custom-made database of mosquito COI genes. Matches were manually queried to the NCBI website to confirm their identity. Mosquito COI sequences are available on NCBI through the GenBank accession numbers KP293419 to KP293425.

### Viral Phylogenetic Analysis and Coverage Estimation

Phylogenetic analysis was performed by comparing the complete and contiguous sequence of each identified viral ORFs to related taxa. Alignments were performed using MAFFT v7.058 and the E-INS-i algorithm (Katoh and Standley, 2013). This alignment algorithm is suitable for sequences that contain multiple conserved regions embedded in long unalignable regions2 . Bayesian analysis was performed using MrBayes v3.1.2 (Ronquist and Huelsenbeck, 2003). The substitution model was determined by allowing MrBayes to sample across the fixed amino acid rate matrices. For each ORF, two independent chains were run for 1,000,000 generations. The resulting average standard deviation of split frequencies is indicated in the caption

2http://mafft.cbrc.jp/alignment/software/algorithms/algorithms.html

to each figure. Tracer v1.5.0 was used to confirm the stationarity of log likelihoods (Rambaut et al., 2013) and the first 25% of the 10,000 total trees were discarded. Results were visualized using FigTree v1.4.0 (Rambaut, 2012; **Figures 2–5**). All trees are available at http://dx.doi.org/10.6084/m9.figshare.1247641. In each case, the alignments contain highly divergent regions without homologs in all taxa, which are encoded with gaps to represent missing data. To confirm these gap sections do not significantly affect reconstruction, we repeated all phylogenetic analyses using only sites that contain data from at least 50% of the taxa. Phylogenetic reconstruction using MrBayes on these shortened alignments found very similar topologies as presented in **Figures 2–5**, and did not change any conclusions (phylogenetic trees resulting from shortened alignments available at http://dx.doi.org/10.6084/m9.figshare. 1247641).

Coverage was determined by mapping the quality checked reads to the nucleotide sequences corresponding to the viral ORFs using Bowtie 2 and default settings (Langmead et al., 2009).

### Results and Discussion

In this study, we used RNA SMS to characterize the microbial communities associated with three mosquito species in northern California: *Culex pipiens*, *Culiseta incidens*, and *Ochlerotatus sierrensis* (previously *Aedes sierrensis*). *Cx. pipiens* (common house mosquito) is a known vector for *West Nile virus* and *St. Louis encephalitis virus* (Farajollahi et al., 2011), *C. incidens* (cool-weather mosquito) is not a major vector for *West Nile virus* or any other disease, and *O. sierrensis* (western treehole mosquito) is a major vector for dog heartworm (Ledesma and Harrington, 2011). Samples were collected from a variety of habitats including sylvatic/wild, rural, suburban, and urban (*sensu* Thongsripong et al., 2013). All three of these species have limited ranges (flight range *<*5 miles), making them good markers for the respective habitats in which they are found. Only female mosquitoes were used, as female mosquitoes, being hematophagous, are the primary vectors of most human and animal diseases.

sequences from this study. This consensus phylogeny was generated in MrBayes (Ronquist and Huelsenbeck, 2003) based on two independent chains run for one million generations to convergence. The average standard deviation of split frequencies was 0.002810. Posterior node probabilities are shown at

uncovered in this study are labeled as such. Closely related mosquito-associated viruses are marked with an asterisk. The accession numbers of all sequences are available at http://dx.doi.org/10.6084/ m9.figshare.1247641.

As many clinically important mosquito-vectored viruses are RNA based (for example, certain *Alphaviruses*, *Bunyaviridae*, and *Flaviviruses*), we focused on the RNA metagenome of these mosquitoes. Nearly 250 million 100 bp paired-end reads were generated on the Illumina HiSeq platform, roughly evenly distributed among the seven libraries (**Table 1**). After quality control and contig assembly, the identity of the viral sequences associated with these samples was determined by querying the translated contigs against a custom database containing the entire NCBI

non-redundant protein database, the PhAnToMe phage protein database (Aziz et al., 2012), the NCBI Viral RefSeq Database, and the *Ae. aegypti*, *An. gambiae*, and *Cx. quinquefasciatus* protein databases (Megy et al., 2012). The identity of the bacterial and fungal sequences associated with these samples was determined by first querying the non-translated contigs to the SILVA SSU and LSU Reference Databases (Quast et al., 2013). Any contigs with a closest match to either bacteria or fungi were then submitted to the RDP Classifier (Wang et al., 2007). Despite using

Ribo-ZeroTM to deplete vertebrate rRNA, adequate bacterial and fungal RNA remained for microbial classification.

included as the outgroup. The average standard deviation of split frequencies

### Viral Communities

Consistent with previous studies (Cook et al., 2013; Chandler et al., 2014; Coffey et al., 2014), viral sequences from the lineages *Bunyaviridae*, *Rhabdoviridae*, and *Narnavirus* were uncovered in these mosquitoes (**Table 2**). In addition, we also detected sequences related to an unassigned double-stranded RNA virus (dsRNA) which has been found to be associated with insects (Spear et al., 2010; Cook et al., 2013). No single viral taxon was found in all samples; however, members of the dsRNA group and the *Rhabdoviridae* were detected in three samples. Interestingly, the San Rafael *Cx. pipiens* sample and the San Francisco *Cx. pipiens* sample each contained representatives of the same three viral taxa – *Bunyaviridae*, *Rhabdoviridae*, and *Narnavirus*. Using the methods employed here, viral sequences were not detected in the Stinson Beach or Bolinas *Cx. pipiens* samples. If viruses are present in these samples, they fall below our detection limit (see Section "Materials and Methods").

sequences are available at http://dx.doi.org/10.6084/m9.figshare.1247641.

### *Bunyaviridae*

In two samples, *Cx. pipiens* from San Francisco and *Cx. pipiens* from San Rafael, we uncovered sequences related to the *Bunyaviridae* L segment. There are five recognized *Bunyaviridae* genera (*Phlebovirus, Hantavirus*, *Nairovirus*, *Orthobunyavirus*, and *Tospovirus*), and all except *Hantavirus* are known to be arthropod-vectored. The *Bunyaviridae* have single-stranded negative-sense genomes consisting of three segments. In the genus *Phlebovirus* (which our sequences are most closely related to, see next paragraph), these segments are designated L (encoding the RNA-dependent RNA polymerase), M (encoding the glycoprotein), and S (encoding the nucleocapsid and an ambisense non-structural gene, the NSs, which is necessary for virulence in vertebrates).

Phylogenetic analysis indicates that the detected *Bunyaviridae* sequences are sister to the *Phlebovirus* genus, which includes the agent of Rift Valley fever, and relatives such as *Gouléako virus* (**Figure 2**). This is in congruence with several other mosquito-associated *Bunyaviridae* that have been found in Thailand (Chandler et al., 2014), France (Cook et al., 2013), and West Africa (Marklewitz et al., 2011). Although the two uncovered sequences are more closely related to each other than to any other taxa, these two sequences differ by approximately 6 amino acids substitutions and a 64 amino acid insertion/deletion at the terminal end (alignments available at http://dx.doi.org/10.6084/m9.figshare.1247641).

Attempts to find the M and S *Bunyaviridae* segments using specific blastp and vFam searches (Skewes-Cox et al., 2014) were unsuccessful (data not shown). While this is surprising given the depth of coverage of the L segment (**Table 2**), we note that a previous metagenomic shotgun sequencing study which found a *Bunyaviridae* L segment likewise did not find either the M or the S segments (Cook et al., 2013). The inability to recover these segments may reflect their greater degree of divergence than the L segment, or their reduced representation in the libraries. Targeted amplification of the missing segments by PCR could compensate for the latter, but may have limited success in the case of divergence. Previous work finds that some mosquito viruses lack the NSs gene, which is one of two genes, along with the nucleoprotein, present on the S segment (Marklewitz et al., 2011; Chandler et al., 2014). The absence of the NSs, which is found in all *Phleboviruses* and is necessary for virulence in vertebrates, and the presence of the nucleoprotein, which is required for replication, is suggestive of an arthropod-only lifecycle in these viruses (Marklewitz et al., 2011; Chandler et al., 2014).

The nucleic acid contigs of both *Bunyaviridae* detected here included stop codons on either side of the ORF suggesting that the entire ORF was found. However, the terminal hairpins that are present in all Bunyaviridae were not recovered, despite implementing PRICE (Ruby et al., 2013) to extend the fragments (data not shown). This suggests that the entire L segment was not present in our data.

### *Rhabdoviridae*

In three of our samples, *O. sierrensis* from Pepperwood Preserve, *Cx. pipiens* from San Rafael, and *Cx. pipiens* from San Francisco, we detected sequences related to the family *Rhabdoviridae*. The *Rhabdoviridae* is a diverse family in the order *Mononegavirales* and consists of six established genera (*Vesiculovirus, Lyssavirus, Ephemerovirus, Novirhabdovirus, Cytorhabdovirus,* and *Nucleorhabdovirus*) along with more than 130 unassigned viruses, such as the *Drosophila*-associated *Sigma viruses* (Kuzmin et al., 2009; Longdon et al., 2010). *Rhabdoviridae* are the causative agents of numerous clinically and economically important diseases of humans, livestock, fish, and plants – some of which are arthropod-vectored (Kuzmin et al., 2009). The genome of the *Rhabdoviridae* consists of a single negative-sense RNA segment with, at minimum, five ORFs corresponding to the nucleoprotein, glycoprotein, phosphoprotein, matrix protein, and the RdRp.

Using our methods, we uncovered partial fragments of the *Rhabdoviridae* nucleoprotein (*Cx. pipiens* from San Rafael and *Cx. pipiens* from San Francisco) and the *Rhabdoviridae* glycoprotein (*O. sierrensis* from Pepperwood Preserve). Since the RdRp is the preferred gene for phylogenetic reconstruction (Bourhy et al., 2005; Longdon et al., 2010), we searched for it within these three libraries using specific blastp and vFam searches (Skewes-Cox et al., 2014), but did not recover the RdRp in our datasets. This is unsurprising since the depth of coverage for the three uncovered fragments was very low compared to many of the other uncovered viruses (**Table 2**).

Phylogenetic analysis was performed using MrBayes on concatenated alignments of the nucleoprotein and glycoprotein fragments against reference sequences for these regions (**Figure 3**). The resulting phylogeny successfully recapitulates the relationships between *Rhabdoviridae* lineages ((((*Vesiculovirus, Ephemerovirus), Sigma viruses*), *Lyssavirus), (Cytorhabdovirus, Nucleorhabdovirus*), *Novirhabdovirus*) found by other studies (Longdon et al., 2010; Coffey et al., 2014). Our mosquitoassociated *Rhabdoviridae* form a clade with *North Creek virus*, isolated from *Cx. sitiens* in Australia (Coffey et al., 2014) and a unnamed virus isolated from *Cx. tritaeniorhynchus* in Japan (Kuwata et al., 2011). This mosquito-associated clade is sister to the *Vesiculovirus, Ephemerovirus,* and *Sigma virus* clades; basal to these is the mosquito-associated *Moussa virus* (Quan et al., 2010).

### Double-Stranded RNA Viruses

In three libraries, we identified sequences that are related to dsRNA viruses previously identified in culicine mosquitoes from France (Cook et al., 2013). These viruses are related to, and have a similar genome structure as, viruses isolated from the plant-feeding insects *Circulifer tenellus* and *Spissistilus festinus* (Spear et al., 2010). The genomes of these viruses consist of a single segment with two ORFs, an RNA-dependent RNA polymerase (RdRp) and a proline–alanine rich protein (PArp) of unknown function. Phylogenetic analysis of the RdRp suggests that the viruses from plant-feeding insects belong to a novel genus of dsRNA viruses (Spear et al., 2010). Phylogenetic analysis using the PArp ORF has not been performed in any previous publication.

We identified sequences related to these dsRNA viruses in three of our samples: *O. sierrensis* from Pepperwood Preserve, *C. incidens* from Mill Valley, and *C. incidens* from San Mateo. Inspired by the discovery of these viral sequences, we searched a prior dataset from our laboratory [*Ae. aegypti* from Thailand (Chandler et al., 2014)] using the methods described above and identified another dsRNA virus. In several instances, a single contig spanned the entire RdRp and PArp ORFs (contig 32 in the *O. sierrensis*sample from Pepperwood Preserve, contigs 84 and 89 in the *C. incidens* sample from Mill Valley, and contig 119 in the *C. incidens* sample from San Mateo; **Table 2**). In the San Mateo sample, one low coverage contig contains partial sequences of both the RdRp and the PArp (contig 1366) and two additional low coverage contigs (contigs 20718 and 26686) contain fragments of the RdRp that partially overlap with, and are nearly identical to, the contig containing the partial sequences.

Phylogenetic analysis was performed using the complete RdRp ORFs and the partial ORF from the San Mateo sample (contig 1366). Viruses from three proposed novel genera, one *Megabirnavirus*, and *Penicillium chrysogenum* virus were included as comparison taxa (Spear et al., 2010; Ito et al., 2013). Our phylogeny reproduces the overall structure of previous work (Ito et al., 2013) and places all mosquito-associated viruses together (**Figure 4**). Interestingly, both the San Mateo and Mill Valley samples have two distinct dsRNA contigs. This suggests that multiple strains of this dsRNA virus are co-circulating within the northern California *C. incidens* population, a situation similar to that occurring in *Circulifer tenellus* in central California (Spear et al., 2013).

In the dsRNA plant-feeding insect viruses, the RdRp and PArp ORFs are not in frame, but are suspected to be transcribed as a fusion via a -1 ribosomal frameshift just prior to the stop codon of the 5 (PArp) ORF. The site of this frameshift is predicted to occur at a G\_GAA\_AAC\_stop motif in the virus of French culicine mosquitoes (Cook et al., 2013). In all of our contigs, the PArp and the RdRp were a single base out of frame. In three of our contigs (84 from Mill Valley *C. incidens*, 32 from Pepperwood Preserve *Cx. pipiens*, and 119 from San Mateo *C. incidens*), we identified a G\_GAA\_AAC\_stop motif and in the remaining contigs (89 from Mill Valley *C. incidens* and 1366 from San Mateo *C. incidens*), we identified a G\_GGA\_AAC\_stop motif. This suggests that the two dsRNA ORFs found here are transcribed as a fusion, although enabled by slightly different motifs.

### *Narnavirus*

Two of our samples, *Cx. pipiens* from San Rafael and *Cx. pipiens* from San Francisco, included sequences that are related to the genus *Narnavirus* in the family *Narnaviridae*. *Narnaviruses* have single-stranded positive-sense genomes consisting of a single ORF that encodes an RdRp (King et al., 2012). The type strain *Saccharomyces 20S virus* infects the Ascomycota *Saccharomyces cerevisiae*, while similar sequences have been found associated with Basidiomycota fungi (Cook et al., 2013) and an oomycete plant pathogen (Cai et al., 2013). While the sister genus to *Narnavirus* is the mitochondrial-infecting *Mitovirus*, the RdRp of plant-infecting *Ourmiavirus* is more closely related to *Narnavirus*, as it was acquired via reassortment (Rastgou et al., 2009), and is therefore used as the outgroup in **Figure 5**.

Phylogenetic analysis finds that the *Cx. pipiens Narnaviruses* are more closely related to the other mosquito-associated *Narnaviruses* (Cook et al., 2013) than to any other taxa (**Figure 5**). Since most of the known *Narnaviruses* are suspected to infect Ascomycete or Basdiocymete fungi [the only exception being the oomycete-infecting *Phytophthora infestans RNA virus* (Cai et al., 2013)], we investigated the presence and identity of potential host fungal sequences in our samples (see Section "Fungal Communities"). We also used the methods described in Section "Sequence Processing and Taxonomy Assignment" to identify fungal sequences in the French mosquito dataset that contained two *Narnavirus* sequences (Cook et al., 2013). We did not find any oomycete SSU or LSU sequences in any of the three datasets (data not shown), thus ruling out that the mosquitoassociated *Narnaviruses* are in fact associated with the presence of oomycetes.

All three mosquito samples positive for *Narnaviruses* also included sequences of Ascomycete fungi. Although no single fungal genus was common to all three mosquito samples, several [such as *Cladosporium* and *Capnodiales* (see Section "Fungal Communities")] were shared between the two northern California samples. Considering that the mosquito-associated *Narnavirus* clade is sister to *Narnaviruses* of fungi, and that fungal sequences were also detected in the mosquitoes, the *Narnaviruses* from mosquitoes reported here and by Cook et al. (2013) may in fact be instances of virus-infected fungal infections of mosquitoes. Alternatively, these mosquito-associated *Narnaviruses* may represent direct infections of mosquitoes. However, *Saccharomyces 20S RNA virus* is vertically transmitted from mother to daughter cells and horizontally through mating (King et al., 2012), suggesting that host switching may be biologically improbable, particularly from a fungal host to a mosquito host. Taken together, our approach may have revealed an interesting scenario: since the method employed here simultaneously characterized both viruses and fungi within a given sample, it can identify patterns of multi-level nested host-parasite associations. In this case, we have identified putative fungal infecting viruses and the fungi they may be infecting, all in the same hosts. This illustrates the potential utility of our broad approach to microbiota characterization – one that could reveal microbial interactions that may ultimately lead to the development of infectious disease control measures.

### Bacteriophages

No phage sequences were uncovered in this dataset, despite including the PhAnToMe phage database in our blastp search. Given that these mosquitoes were associated with a diversity of bacteria (see Section "Bacterial Communities") and that the microbiota of other animals consists of abundant phage (De Paepe et al., 2014), it is reasonable to presume that phages are indeed associated with these mosquitoes as well. There exist two explanations for our inability to find phage associated with these samples. First, our dataset was RNA based and most phages are DNA based. Indeed, known phages of the two most widespread bacteria associated with these samples (*Bacillus* and *Escherichia/Shigella*; see Section "Bacterial Communities") are double-stranded DNA based. Second, perhaps the phages infecting the bacteria within these mosquitoes are undescribed and therefore would not be represented in the databases.

### Bacterial Communities

Contigs that could be confidently identified as bacteria using the RDP Classifier were found in all seven samples (**Table 3**). Since there are known limitations of using rRNA as a measurement of either microbial abundance or activity (Blazewicz et al., 2013), the quality-checked reads were not mapped back to these contigs to determine the number of reads used to build each contig. Therefore the data presented here (and of the fungi discussed below) should only be interpreted as representing the presence of that particular taxon and not that taxon's abundance within the community.

Contigs (those that passed the initial 300 bp cutoff) ranged from 328 to 1546 bases in length, which means that some cover nearly the entire 16S gene and all cover at least one hyper-variable region. A single contig with a bacterial match (*Bacillus* in the Mill Valley *C. incidens* sample) was found to be chimeric (data not shown) and only the bacterial portion was used for classification. We identified the following bacterial genera in our samples: *Actinobacillus, Anaplasma* [*Wolbachia* (see below)]*, Bacillus, Corynebacterium, Escherichia/Shigella, Haemophilus, Moraxella, Propionibacterium, Pseudomonas* and *Wandonia*. When a genus was found multiple times in the same sample (for example, *Moraxella* in the Mill Valley *C. incidens* sample), the contigs spanned different regions of the 16S gene (see blastn results at http://dx.doi.org/10.6084/m9.figshare.1247641).

*Bacillus* and *Escherichia/Shigella* were the most widespread genera, as they were found in all seven samples. Notably, a recent review by Minard et al. (2013) concludes that female mosquitoes are mostly colonized by Gammaproteobacteria (of which *Escherichia/Shigella* is a member) and that males are dominated by Firmicutes, such as *Bacillus*. The results of our study, in which only females were investigated, are therefore in partial agreement with this conclusion.

*Anaplasma* is a genus in the order Rickettsiales. *Wolbachia* is a notable and widespread arthropod-associated bacteria in this order (Hilgenboecker et al., 2008) and we therefore investigated the *Anaplasma* contigs further. Querying these contigs to the NCBI database finds them to be greater than 98% identical to *Wolbachia pipiens* (data not shown). The two samples containing *Wolbachia* were *Cx. pipiens* from San Rafael and Stinson beach. The third *Cx. pipiens* sample, from Bolinas, also had contigs identified as *Anaplasma*, although their confidence was below the 90% cutoff in RDP (**Table 3**). Querying these contigs to NCBI found they were greater than 99% identical to *Wolbachia pipiens* (data not shown, all contigs are available at http://dx.doi.org/10.6084/m9.figshare.1247641). Inspection of the contigs from *C. incidens* and *O. sierrensis* that were shorter than 300 bp or below the 90% confidence cutoff did not reveal any *Anaplasma* or *Wolbachia*-like sequences.

*Wolbachia* was first described in *Cx. pipiens* mosquitoes (Hertig and Wolbach, 1924). A survey of California mosquitoes (including *Culex, C. incidens,* and *Ochlerotatus*) for *Wolbachia* infections found that, out of 296 individuals and 14 species, only *Cx. pipiens* mosquitoes were infected (Rasgon and Scott, 2004). How *Wolbachia*, which is vertically transmitted, spreads into new taxa is not well understood, but its absence in certain species could be due to inhibition by other members of the mosquito microbiota, as has been shown in *Anopheles* mosquitoes (Hughes et al., 2014). While this study did not find any consistent microbiota differences between *Wolbachia*-infected and *Wolbachia*-free mosquitoes, more comprehensive studies may find such differences that might help account for natural *Wolbachia* distributions.

Of final note, two of the bacterial genera identified, *Propionibacterium* and *Corynebacterium,* are commonly associated with human skin (Grice and Segre, 2011). There are two potential explanations for their presence in these samples. One is that they were acquired from human hosts during a blood meal. Alternatively, these could be contaminants introduced during collection or processing for sequencing. Future studies will include a negative control to identify this possibility (Salter et al., 2014).

### Fungal Communities

Contigs that could be confidently classified as fungi using the RDP Classifier were found in all samples except *O. sierrensis* from Pepperwood Preserve (**Table 4**). Manually querying the putative fungal contigs (identified by the initial SILVA database search) from this sample against the NCBI database confidently identified a fungus that was nonetheless below the RDP confidence cutoff (**Table 4**).

Fungal contigs ranged in length from 302 to 1173 bp. Fungi from 13 families were present in our data (**Table 4**). In contrast to the bacterial data, there is no single fungal genus that is present in all samples and only two fungi were found in multiple samples. Specifically, *Chromocleista* was found in *Cx. pipiens* from both Stinson Beach and Bolinas, and *Cladosporium* was found in *Cx. pipiens* from both San Rafael and San Francisco. All the remaining fungal taxa identified were unique to a single sample. More in-depth sequencing may reveal widespread fungi that exist at numbers within the host too low for our methods to detect.

Some of the fungi identified within these samples may not exist in an intimate relationship with mosquitoes. For example, several taxa within the Basidiomycota have a characteristic multicellular stage, such as the *Mycenella* gilled mushrooms identified in the *Cx. pipiens* sample from Stinson Beach. Most likely, dispersing spores of these fungi contaminated water sources used by adult mosquitoes as breeding sites. Similarly, the fungal taxa *Malassezia*, identified in *Cx. pipiens* from San Francisco, is a common human commensal and opportunistic pathogen (Gupta et al., 2004) that likely contaminated the sample through human contact, as with the bacteria discussed above.

Much previous work has focused on using entomopathogenic fungal taxa for mosquito control (Scholte et al., 2004; de Faria and Wraight, 2007). Interestingly, several of the fungi identified belong to families that include many known entomopathogenic fungi including Clavicipitaceae (found in the *Cx. pipiens* sample from San Francisco) and Cordycipitaceae (found in the *O. sierrensis* sample from Pepperwood Preserve). For future mosquito collections, histologic examination would be necessary to identify visible signs of disease in positive mosquitoes.

This study represents, to our knowledge, the first to use culture-independent techniques to characterize fungal communities in natural populations of mosquitoes. As with other animal hosts [e.g., humans (Huffnagle and Noverr, 2013) and *Drosophila* (Broderick and Lemaitre, 2012)], fungi associated with mosquitoes are relatively understudied compared to bacteria. One study of laboratory-raised *Ae. aegypti* identified *Saccharomyces, Penicillium,* and *Aspergillus* using methods similar to those employed here (Bishop-Lilly et al., 2010). While the current study did find these, or closely related, taxa, more work is needed to fully determine how the fungal communities of laboratory-raised and wild-caught mosquitoes compare.

### Mosquito Identification

Five of the seven mosquitoes were identified morphologically or by PCR amplification and Sanger sequencing of the COI gene prior to SMS. For the two specimens that were not identified beforehand, contigs matching mosquito COI and mosquito LSU and SSU were queried against the NCBI database and indicated that these mosquitoes were *Cx. pipiens* and *O. sierrensis* (data not shown, all contigs available at http://dx.doi.org/10.6084/m9.figshare.1247641).

### Conclusion

In this study, we used SMS to identify viral, bacterial, and fungal sequences associated with naturally collected mosquitoes. Notably, in all three cases of mosquito-associated viruses (i.e., the *Bunyaviridae*, *Rhabdoviridae*, and dsRNA viruses) the viruses described here were more closely related to other mosquitoassociated viruses than to any other viruses. This is despite the fact that the mosquitoes belonged to a variety of taxa (e.g., *Ae. aegypti, C. incidens, O. sierrensis,* and several species of *Culex*) and were collected from a variety of locations (e.g., northern California, Thailand, France, western Africa, and Australia). Furthermore, it should be noted that two of the viral groups discussed here (*Bunyaviridae* and *Rhabdoviridae*) have been extensively studied due to their importance in human and animal health and are therefore well characterized in terms of host breadth. Taken together, a pattern of family level host-specificity (i.e., to the family Culicidae) is emerging. We look forward to future research to sample the viral taxa more extensively and help answer the question of host-specificity and host-switching in the viral tree of life.

### References

With respect to the bacterial and fungal components of the mosquito microbiota, two opposing patterns emerged: several bacterial taxa were widespread regardless of mosquito species or habitat, whereas fungal communities were distinct between individuals even within host species or habitat. These patterns may be an artifact of insufficient sample size and sequencing depth. Much of our sequencing effort was needlessly expended on host rRNA, as over a quarter of the total data can be attributed to the mosquito ribosomal large subunit or small subunit. However, recent techniques have been developed to selectively remove sample-specific rRNA from environmental (Stewart et al., 2010) and mosquito (Kukutla et al., 2013) samples. The use of these techniques will greatly increase the power to detect bacteria, fungi, and viruses in field-caught mosquitoes. Furthermore, many more samples could be combined into a single experiment to allow for the inclusion of multiple replicates per population, host species, or habitat type. Additionally, the physical environment of the mosquitoes, such as larval breeding sites and attractants used in collection, can be sampled. In this way, the complex interactions between host, environment, and all the microbes inhabiting both, can be determined.

### Acknowledgments

We thank the California Academy of Sciences' Student Science Fellows high school enrichment program for support, especially Roberta Brett, Sean Edgerton, and Ryan Hulett. This program was generously supported by the Gordon and Betty Moore Foundation, with additional support provided by Mary Austin and Brewster Kahle. V. Piper Kimball and team members of the Marin/Sonoma Mosquito and Vector Control District executed mosquito collections at Stinson Beach, Bolinas, and San Rafael. Durrell Kapan assisted in other mosquito collections. Panpim Thongsripong, Barbara Arrighi, and Anna Sellas provided technical assistance. We thank Jonathan Eisen, Jack Dumbacher, and Joseph DeRisi for use of their servers.

### Supplementary Material

The supplemental data associated with this work is available at: http://dx.doi.org/10.6084/m9.figshare.1247641. This includes alignments, tree files, blast outputs, and other intermediate files.

clone libraries contain chimeras. *Appl. Environ. Microbiol.* 72, 5734–5741. doi: 10.1128/AEM.00556-06

high throughput sequencing. *PLoS ONE* 6:e24758. doi: 10.1371/journal.pone. 0024758

viruses and seasonal prevalence of known viruses, Nosema, and Crithidia. *PLoS ONE* 6:e20656. doi: 10.1371/journal.pone.0020656

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Copyright © 2015 Chandler, Liu and Bennett. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

## Metagenomics of plant and fungal viruses reveals an abundance of persistent lifestyles

### *Marilyn J. Roossinck\**

*Department of Plant Pathology and Environmental Microbiology, Center for Infectious Disease Dynamics, Pennsylvania State University, University Park, PA, USA \*Correspondence<EMAIL_ADDRESS>
#### *Edited by:*

*Bas E. Dutilh, Radboud University Medical Center, Netherlands*

#### *Reviewed by:*

*Hano Maree, Agricultural Research Council, South Africa Huiquan Liu, Northwest A&F University, China*

**Keywords: virus ecology, mycoviruses, dsRNA viruses, persistent plant viruses, mutualistic viruses**

### **INTRODUCTION**

Most studies of plant viruses have focused on the acute viruses that cause disease in crop and ornamental plants. These viruses are transmitted horizontally, often by insect vectors, and are occasionally transmitted vertically. Although known for at least four decades, the persistent viruses of plants are very poorly studied. These viruses were previously called "cryptic" because they did not appear to illicit any symptoms in infected plants (Boccardo et al., 1987). Persistent plant viruses are not known to be transmitted horizontally, although phylogenetic evidence suggests some level of transmission (Roossinck, 2010). They are vertically transmitted at nearly 100% levels through both ova and pollen (Valverde and Gutierrez, 2007). They have been identified in metagenomic studies by their similarity to known persistent viruses, and because they lack any movement protein, a feature of all known acute viruses that must move through the plant plasmodesmata to establish a systemic infection. Persistent viruses do not move between plant cells, but rather infect every cell and move by cell division.

Most plant persistent viruses have double-stranded (ds) RNA genomes, and encode only an RNA dependent RNA polymerase (RdRp) and a coat protein. Of the well-characterized persistent plant viruses, those in the *Endornaviridae* are the exception. These viruses have a singlestranded (ss) RNA genome, based on their RdRp, and encode a large polyprotein that does not have any apparent coat protein, but encodes a number of additional domains that appear to be derived from diverse sources (Roossinck et al., 2011). They are usually found as dsRNA replicative intermediates.

Viruses of fungi have very similar lifestyles to plant persistent viruses, and several virus families are shared between plants and fungi. Phylogenetic evidence indicates that virus transmission has occurred within and between the two kingdoms (Roossinck, 2010; Roossinck et al., 2011).

Fungal viruses are even less wellstudied than plant viruses, and the diversity of these viruses remains mostly unknown. A majority of known fungal viruses have dsRNA genomes, some have ssRNA genomes, and a few examples of DNA viruses are known (Yu et al., 2010). Recently a negative sense ssRNA virus was characterized from a fungus (Liu et al., 2014). Similar to plant viruses, most fungal viruses have been studied in the context of pathogenic fungi. The discovery of the hypovirulence phenotype of *Cryphonectria hypovirus* 1 that suppresses the disease phenotype of the chestnut blight fungus led to a search for other examples that could be exploited to mitigate the effects of plant pathogenic fungi [reviewed in Dawe and Nuss (2013)].

### **VIRUS DISCOVERY IN PLANTS AND FUNGI**

Deep sequencing is proving to be a useful technique for just about everything these days, and the methods have been applied to metagenomic studies of viruses. Unlike studies of other microbes, viruses cannot be analyzed through the use of any universal conserved sequences or motifs, and a variety of techniques have been employed to enrich for viral nucleic acids before sequence analysis. Studies in aquatic viruses have been reported for a number of years (Angly et al., 2006; Labonté and Suttle, 2013). More recently plant viruses have been studied through metagenomics as well (Roossinck, 2012; Stobbe and Roossinck, 2014). A large variety of studies have been done on many different scales, from individual plants to ecosystems. In some studies a single plant species has been targeted, in others a broader sweep is used. These studies are explored in detail in a review by this author and others to be published elsewhere. Here I will explore the discovery of persistent viruses that are extremely common in plants and fungi, but poorly studied, and discuss the implications of these viruses in the ecology of plants and fungi.

Different methods of detection have yielded different levels of persistent viruses. Use of dsRNA-enriched samples yielded very high levels of persistent viruses in plants (Roossinck, 2012). Using the small RNAs involved in plant immunity (siRNAs) has been less successful at detecting many persistent viruses in plants. While the complete sequence of a known endornavirus was assembled with siRNAs (Sela et al., 2012), no novel endornaviruses have been reported with this method. A few siRNAs have been found for partitiviruses, chrysoviruses, and totiviruses, but with very limited genome coverage (Kreuze, 2014). It is likely that these viruses are not subjected to silencing; with the exception of endornaviruses, they do not expose their dsRNA to the cell, but rather retain their genomes within the virions and extrude only ssRNA into the cytoplasm (Safari and Roossinck, 2014).

Virus discovery in fungi is very limited. Most analyses have been done on fungi of economic importance such as plant pathogenic fungi. A survey of viruses from endophytic fungi derived from two plant species in a wild plant community indicated that the diversity of viruses in this system was much greater than the diversity of fungi, which was in turn much greater than the diversity of plants (Feldman et al., 2012).

In most cases of fungal virus studies, viruses have been discovered from cultured fungi. This eliminates the majority of fungi, which are not culturable (Blackwell, 2011), but have been discovered from environmental samples through specific gene analysis such as ribosomal RNA-related regions and other genes (Seifert, 2009). Traditionally fungi acquired from nature are "purified" by single spore isolation. These practices result in a gross under-estimate of fungal viruses, as many viruses are lost during culture, especially on solid media (unpublished observation), and single spore isolation is a common strategy to obtain cultures "cured" of their viruses. Even though next generation sequencing methods allow for deeper analysis of environmental samples, finding new viruses in fungi without culture is technically challenging. Although some reports have indicated that fungal viruses can be shed into the media when cultured, there is little evidence of extracellular accumulation of most fungal viruses. The lack of any conserved sequences in viruses, such as house-keeping or barcoding genes found in all other life forms, means that sequence-specific primers cannot be used. For viruses of endophytic fungi that have been the focus of fungal virus research in the author's lab, the minimal amount of fungal tissue in plants makes any analysis nearly impossible without culturing the fungus out of the plant. Hence for now we must settle for this very low estimate of fungal virus diversity.

### **COMMON THEMES FROM PLANTS AND FUNGI**

Plants and fungi share several families of viruses. The *Partitviridae* and the *Endornaviridae* are recognized by the International Committee for the Taxonomy of Viruses (King et al., 2012) as infecting both plants and fungi, but biodiversity surveys of plant viruses have also identified members of the *Totiviridae* and *Chrysoviridae* families that traditionally are considered fungal viruses (Roossinck, 2012), and a chrysovirus was recently characterized from radish (Li et al., 2013). In fact viruses from these and related families make up over half of the viruses identified in wild plants (Roossinck, 2012). In plants these viruses appear to maintain a persistent lifestyle (Roossinck, 2010), remaining associated with their hosts for many generations with nearly 100% vertical transmission. Less is known about the lifestyles of fungal viruses. There are few reports of truly acute viruses in fungi. Recently a DNA virus from *Sclerotinia sclerotiorum* was shown to be infectious as a purified virus particle, although it is not clear if this is a mechanism for transmission in nature (Yu et al., 2013). Unlike the plant persistent viruses, fungal viruses can be transmitted between closely related strains of fungus through anastomosis (Milgroom and Hillman, 2011), and evidence of cross-species transmission is apparent in phylogenetic analyses of *Cryphonectria hypovirus* (Liu et al., 2003) and partitiviruses in the *Heterobasdion* (Vainio et al., 2011).

Persistence and high levels of vertical transmission in parasites are correlated with commensal or mutualistic lifestyles (Villarreal, 2007; Márquez and Roossinck, 2012). In some cases we know that persistent viruses are mutualistic (Nakatsukasa-Akune et al., 2005; Márquez et al., 2007). In many cases we don't know enough about their biology to assess their symbiotic lifestyle, but in plants few have any evidence of negative effects on their hosts. For two persistent viruses in *Heterobasidion* species, virus lifestyle was dependent on the fungal environment (Hyder et al., 2013). A complicating factor in understanding the ecology of persistent viruses is that the well-studied persistent viruses are mainly from crop plants, or from economically important fungi; there is virtually no information about any roles these viruses may play in the natural environment where the virus-host relationships evolved.

### **CONCLUSIONS**

The abundance of persistent viruses in plants and fungi imply functions that may contribute to the biology of the host. Unfortunately we have little ecological data about these viruses, and since they often cause no disease they have not been the subject of intensive study. Datamining from transcriptomic, genomic and metagenomic studies may allow us to address the true ecological role of these viruses. For example, the partitiviruses have poly-A tails, and may be detectable in transcriptome analyses (Jiang et al., 2013). In some cases persistent virus sequences are found integrated into plant or fungal genomes (Liu et al., 2010; Chiba et al., 2011). Deeper analyses along these lines may provide data on time-lines of persistent virus-host relationships.

### **ACKNOWLEDGMENTS**

The author acknowledges the National Science Foundation grant numbers EF-0627108, EPS-0447262, IOS-0950579, and IOS-1157148 and The United States Department of Agriculture grant number OKLR-2007–01012 for previous and continuing research support, and the Pennsylvania State University.

### **REFERENCES**

**Conflict of Interest Statement:** The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 04 December 2014; accepted: 16 December 2014; published online: 12 January 2015.*

*Citation: Roossinck MJ (2015) Metagenomics of plant and fungal viruses reveals an abundance of persistent lifestyles. Front. Microbiol. 5:767. doi: 10.3389/fmicb. 2014.00767*

*This article was submitted to Virology, a section of the journal Frontiers in Microbiology.*

*Copyright © 2015 Roossinck. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# **Tales from the crypt and coral reef: the successes and challenges of identifying new herpesviruses using metagenomics**

### *Charlotte J. Houldcroft 1\* and Judith Breuer 1,2*

*1 Infection, Inflammation and Rheumatology, Institute of Child Health, University College London, London, UK, 2 Division of Infection and Immunity, University College London, London, UK*

#### *Edited by:*

*Richard J. Hall, Institute of Environmental Science and Research, New Zealand*

#### *Reviewed by:*

*Marta Canuti, Memorial University of Newfoundland, Canada Jenny L. Draper, Ministry for Primary Industries, New Zealand*

#### *\*Correspondence:*

*Charlotte J. Houldcroft, Infection, Inflammation and Rheumatology, Institute of Child Health, University College London, Cruciform Building, 30 Guilford Street, London WC1E 6BT, UK<EMAIL_ADDRESS>
#### *Specialty section:*

*This article was submitted to Virology, a section of the journal Frontiers in Microbiology*

> *Received: 03 December 2014 Paper pending published: 11 January 2015 Accepted: 20 February 2015 Published: 13 March 2015*

#### *Citation:*

*Houldcroft CJ and Breuer J (2015) Tales from the crypt and coral reef: the successes and challenges of identifying new herpesviruses using metagenomics. Front. Microbiol. 6:188. doi: 10.3389/fmicb.2015.00188* Herpesviruses are ubiquitous double-stranded DNA viruses infecting many animals, with the capacity to cause disease in both immunocompetent and immunocompromised hosts. Different herpesviruses have different cell tropisms, and have been detected in a diverse range of tissues and sample types. Metagenomics—encompassing viromics—analyses the nucleic acid of a tissue or other sample in an unbiased manner, making few or no prior assumptions about which viruses may be present in a sample. This approach has successfully discovered a number of novel herpesviruses. Furthermore, metagenomic analysis can identify herpesviruses with high degrees of sequence divergence from known herpesviruses and does not rely upon culturing large quantities of viral material. Metagenomics has had success in two areas of herpesvirus sequencing: firstly, the discovery of novel exogenous and endogenous herpesviruses in primates, bats and cnidarians; and secondly, in characterizing large areas of the genomes of herpesviruses previously only known from small fragments, revealing unexpected diversity. This review will discuss the successes and challenges of using metagenomics to identify novel herpesviruses, and future directions within the field.

#### **Keywords: metagenomics, virus discovery, next-generation sequencing, herpesviruses,** *de novo* **assembly**

### **Introduction**

Herpesviruses are double-stranded DNA viruses which infect many members of the animal kingdom, including mammals, birds and reptiles (McGeoch and Gatherer, 2005), fish (Lepa and Siwicki, 2012), amphibians (Davison et al., 2006), molluscs (Nicolas et al., 1992; Savin et al., 2010), and cnidarians (Vega Thurber et al., 2008; Grasis et al., 2014). They typically infect a large proportion of their target population, spreading through a variety of horizontal and vertical routes. They can cause disease in many settings, e.g., hemorrhage disease in elephants caused by elephant endothelial herpesviruses (Wilkie et al., 2014); cancers caused by the human gammaherpesviruses Epstein-Barr virus (EBV) and Kaposi's sarcoma-associated herpesvirus (KSHV) (Taylor and Blackbourn, 2011); pulmonary infection in cats caused by feline herpesvirus 1 (Thiry et al., 2009) and terrapene herpesvirus 1-associated pneumonia in Eastern box turtles (Sim et al., 2014); seasonal mass mortality in oysters (Nicolas et al., 1992); and herpes simplex virus (HSV) encephalitis in humans (Whitley, 2006). Disease may be associated with primary infection, reactivation of latent infection, immune

**Abbreviations:** EEHV, elephant endotheliotropic herpesviruses; LCV, lymphocryptovirus.

suppression, or immune senescence. Previous approaches to herpesvirus discovery have utilized a range of methods, discussed in greater depth elsewhere (Bexfield and Kellam, 2011).  Metagenomics is the newest approach to add to this list: a deep-sequencing-based analysis of the nucleic acid contents of a cellular or acellular sample, unreliant on tissue culture and without reference to prior knowledge of which viruses are present in a sample.

Because herpesviruses are widespread among animals, herpesvirus-like sequences are likely to be present in many metagenomic studies. There are also features of herpesvirus molecular biology which increase the likelihood that deepsequencing studies which are not explicitly metagenomic in nature will nevertheless detect herpesviruses, which may or may not be disease-associated. As widespread DNA viruses, any animal genome sequencing project that uses DNA from primary tissue or samples (e.g., saliva or blood) is likely to include sequences from herpesviruses present in that organism. Even though virus detection is not a primary target of host genome sequencing studies, viruses present in the sample will form a proportion of the sequence reads, and herpesviruses have been discovered in exactly such data (Aswad and Katzourakis, 2014). In this review article, we summarize some of the success stories in identifying novel herpesviruses through metagenomics, and offer a warning on the topic of virus discovery.

## **Finding the Needle in a Haystack**

Much herpesvirus genome sequencing in the past has relied upon large volumes of DNA, acquired through culturing the virus in permissive cell lines (Baer et al., 1984; Lei et al., 2013; Lin et al., 2013). In uncultured or unenriched samples, herpesvirus DNA is present in much smaller proportions than host DNA (Depledge et al., 2011). Herpesvirus genomes are typically hundreds of kilobases in length, complicating the task of Sanger or ultra-deep sequencing of herpesvirus genomes. The problem of poor virus to host DNA ratios applies equally to herpesvirus discovery. It is perhaps unsurprising that many herpesviruses are identified using consensus or degenerate PCR primer sets, but their genomes remain unsequenced (e.g., Bodewes et al., 2014; Sim et al., 2014). The success of identifying novel herpesviruses or of characterizing their genomes further is influenced by the similarity of the novel herpesvirus to related viruses. Some uncharacterized herpesviruses may map well to related herpesviruses, and are relatively easy to assemble from deep-sequence data, or design consensus primers to amplify. Divergent herpesviruses are more difficult to detect and assemble using traditional methods, mapping poorly to existing reference sequences (Pop, 2009). Highly similar and divergent herpesviruses alike may require *de novo* assembly to generate accurate consensus sequences, and poten-

## **Deep-Sequencing of Deltaherpesviruses**

Elephant endotheliotropic herpesviruses 1A (EEHV1A) and 1B (EEHV1B), causes of hemorrhage disease in African and Asian elephants, required deep sequencing and *de novo* assembly to characterize their complete genome sequences (Wilkie et al., 2013). EEHVs cannot yet be propagated in tissue culture and must be sequenced from primary tissue—in this case, necropsy tissue from a juvenile Asian elephant.While EEHVs were known of from traditional PCR and Sanger-based methods, the complete genome sequence proved to be divergent from known herpesviruses, with 60 novel herpesvirus genes identified in the EEHV1A genome (Ling et al., 2013). Whole genome sequencing of EEHV5 revealed further genetic diversity between strains, with only 60% similarity between EEHV1A and EEHV5 (Wilkie et al., 2014), and 25% nucleotide divergence between coding regions of EEHV1A and EEHV2 (Richman et al., 2014). The use of deep-sequencing to resolve the genome sequences of EEHVs expands our knowledge of herpesvirus diversity, and may lead to the identification of putative deltaherpesviruses in other species.

## **Herpesviruses Hiding in Plain Sight**

Public sequence repositories containing mapped and unmapped reads from large genome-sequencing projects are fertile territory for identifying novel herpesviruses. Untargeted deep and ultradeep sequencing of mammalian genomes from primary samples can serve as metagenomic datasets for researchers able to reanalyze the data in ways that may not have been anticipated in the original project.Aswad and Katzourakis (2014)screened 14 whole genomes from nine primate sequences available in GenBank for known herpesvirus protein sequence clusters, hoping to identify novel herpesviruses. Characterizing primate herpesvirus diversity is of interest to researchers for a number of reasons: herpesviruses have both zoonotic (Hummeler et al., 1959) and anthroponotic (Huemer et al., 2002) potential in primates, with high fatality rates for humans and non-human primates (Estep et al., 2010); increasing our knowledge of herpesviruses in other species may provide new model systems for understanding how herpesviruses cause disease in humans (Staheli et al., 2014); and the identification of putative endogenous herpesviral elements within primate genomes serves as a fossil record for the herpesviridae (Aswad and Katzourakis, 2014).

Aswad and Katzourakis successfully identified novel herpesvirus-like sequences, and constructed partial herpesvirus genomes, from the whole genome sequences of two primates not previously known to carry herpesviruses: the aye-aye (*Daubentonia madagascariensis*) and Philippine tarsier (*Tarsius syrichta*). Aswad and Katzourakis analyzed both mapped and unmapped reads from the primate genomes, enabling them to identify both exogenous, and chromosomally integrated, potentially endogenised, herpesviruses. The tarsier herpesvirus sequences were particularly interesting because they are the first report of a potentially endogenous herpesvirus, related to human herpesvirus 6B which chromosomally integrates in a small proportion of humans (Luppi et al., 1993; Tanaka-Taya et al., 2004). The tarsier herpesvirus sequences they reported were heavily mutated, not replication competent and contiguous with stretches of host DNA, supporting the view that this virus has become fully endogenised in tarsier evolutionary history. They were also able to assemble much larger regions of the bonobo (*Pan paniscus*) herpesvirus *Pan paniscus lymphocryptovirus 1* (PpanLCV1), covering an estimated 45% of its genome (*∼*78,000 bp), which was only previously known from small fragments (*∼*3000 bp; Ehlers et al., 2003).

While a successful strategy for discovering novel herpesviruses, this approach differs from that employed by metagenomics studies. The authors acted with a reasonable assumption that herpesviruses might be present in the data they were screening, searching only for sequences with similarity to previously reported herpesvirus-like sequences.

Successful identification of novel herpesviruses in host genome sequence datasets also relies on the genomes being sequenced from primary tissue (blood and liver biopsy, in two cases) rather than cell lines, which do not contain the diversity of viruses found in a primary tissue sample. For example, the bonobo genome was sequenced from a bonobo lymphoblastoid cell line which was immortalized in the laboratory using the human herpesvirus EBV (Prufer et al., 2012); this dataset is therefore much more likely to contain human EBV sequences, present in every cell, than primate gammaherpesviruses. Utilizing uncultured primary samples is at the heart of metagenomics.

### **Going Batty for Metagenomics**

Metagenomics may be seen as the opposite side of the coin to host genome sequencing projects: while host genome sequence assembly discards all sequences which are not on target to the host genome, metagenomic studies focus on non-host sequences (reviewed in Mokili et al., 2012). Eliminating or reducing the amount of host nucleic acid present in sequencing data sets has been addressed in some metagenomics studies of bat viral diversity, with success—three recent studies of chiropteran (bat) metagenomes have identified novel herpesvirus sequences, findings confirmed by a fourth study. Bats represent a fifth of the classified mammalian species on earth, and are thought to be a significant reservoir of emerging infectious diseases, particularly those bats found in urban areas or used by humans as a food source (Luis et al., 2013).

Straw-colored fruit bats (*Eidolon helvum*) are widely distributed across Africa and eaten as bushmeat—and therefore a potential zoonotic reservoir (Baker et al., 2013). Studying the viral diversity of this bat is therefore of significant interest to researchers. The sampling approach taken by Baker et al. (2013) was to pool the nucleic acids taken from three sample types: urine, collected non-invasively, and throat and lung swabs obtained from wild, anaesthetized bats from rural and urban locations in Ghana, with sequences assembled *de novo*. They identified sequences from previously unknown bat herpesviruses, with most reads belonging to beta and gammaherpesviruses. The largest proportion of herpesviral reads were obtained from throat swabs, underlining the importance of choosing appropriate primary samples in herpesvirus discovery. The presence of novel herpesviruses was then confirmed by PCR amplification of fragments from the original samples.

Multiple novel herpesviruses were also discovered in a sample of bats from China, which included a broader sampling strategy than in the Baker et al. (2013) study. Anal and pharyngeal swaps were collected from 11 bat species (Wu et al., 2012). In order to maximize the ratio of viral to host/eukaryotic nucleic acid, samples were passed through a filter to remove bacterial and eukaryotic cells. Sequence-independent nucleic acid amplification was used to increase the available genetic material for sequencing, and deep-sequence reads were aligned to the NCBI nt database, identifying novel bat herpesviruses. Even after physical filtering, only 1.2% of sequencing reads were viral, but among these reads, they were able to identify sequences from two previously unknown betaherpesviruses and two gammaherpesviruses. Wu et al. (2012) then used PCR amplification and genome walking to confirm the presence of novel herpesviruses and obtain longer sequences for phylogenetic analysis from the original sample. Clearly low herpesvirus abundance in the initial sample was a factor in preventing whole-genome assembly from the primary material; but metagenomics was integral to identifying the viruses.

A similar study of seven bat species from eastern North American detected a further novel bat betaherpesvirus (Donaldson et al., 2010), from saliva and fresh feces. Samples were passed through a filter to reduce non-viral nucleic acid carry-over, and pooled by species, sex, and age for sequencing to increase the total pool of nucleic acid available. Sequence-independent sample amplification with random hexamers was used, followed by deep-sequencing and *de novo* assembly.

A metagenomic study of French bats (Dacheux et al., 2014) which examined lung, liver and brain tissue from the carcasses of five species was also able to find evidence of the herpesviruses reported by previous studies, reinforcing the ubiquity of herpesviruses among the chiropterans. It is further evident from these studies that saliva, throat swabs and primary tissue samples are integral to herpesvirus discovery.

### **An Ocean of Herpesviruses**

Metagenomics is changing our ideas of where future herpesviruses may be discovered, providing evidence of herpesvirus-like DNA in cnidarians. A metagenomic study of a salt-water coral virome revealed herpesvirus-like sequences in *Porites compressa*. The researchers performed deep-sequencing of the coral in optimal and "stressed" conditions, including increased temperatures, increased acidity, and nutritional stress. They found a small number of herpesvirus-like reads in the coral harvested under optimal conditions, but a much greater number of herpesvirus-like reads were identified in the stressed coral samples. The researchers suggest this is due to lytic reactivation of herpesviruses within the coral under conditions of stress, reflecting what is known about herpesvirus reactivation in mammals. They were able to PCR amplify and sequence a gene with moderate identity to the thymidylate synthase gene from Herpesvirus saimiri 2 (Vega Thurber et al., 2008).

The same research group followed this up with a metagenomic analysis of four species of coral (*Acropora*, *Diploria*, *Montastraea*, and *Porites*) to test for associations between herpesviruslike sequence abundance and disease in cnidarians (Soffer et al., 2014), and used transmission electron microscopy to identify herpesvirus-like particles within the cells of healthy coral. Using *de novo* assembly, they found that herpesvirus-like sequences were more common in healthy coral than diseased coral. This result contradicts the previous study, and highlights the difficulties in making quantitative comparisons between metagenomic datasets, when sequences from novel viruses cannot be normalized to gene or genome length for quantification in the manner of RNA-seq data (Mortazavi et al., 2008). While herpesvirus-like sequences had previously been detected in molluscs (e.g., Davison et al., 2005), metagenomics was instrumental in finding herpesviruslike sequences in cnidarians. Subsequent metagenomic studies of fresh water cnidarians have also found herpesvirus-like sequences in a number of species of hydra (Grasis et al., 2014).

### **Novel Virus or Novel Contaminant?**

There are warnings attached to highly sensitive deep-sequencing and metagenomic discovery of any virus: principally, the risk of contamination. Deep-sequencing and *de novo* assembly of viral nucleic acids provide unprecedented opportunities to identify new viruses, some of which may be pathogenic, but these technologies may also detect contaminating nucleic acids that can be introduced to the sample extraction and sequencing pipeline at many different points.

There have been a number of high profile reports over the last 5 years of "novel" viruses (Xu et al., 2013), which have later been identified as nucleic acid contaminants from commercial sequencing products (e.g., Naccache et al., 2013, 2014; Rosseel et al., 2014). The origins of these nucleic acids are as diverse as mice (Erlwein et al., 2011) and algae (Naccache et al., 2013, 2014). Viral contamination of cell lines (Hue et al., 2010), especially those which have been xenografted during their life in tissue culture (Griffiths et al., 1997, 2002; Paprotka et al., 2011), has also proven to be a significant problem. The story of these "rumor viruses" (Weiss, 2010) must serve as a constant reminder to the virus discovery and metagenomics community. Furthermore, the problem is not confined to viruses, encompassing bacterial contamination of nucleic acid extraction spin columns (Salter et al., 2014) and other laboratory reagents. Screening of control samples and potentially of common laboratory reagents used as each point of the sample preparation process may be necessary for researchers to be confident that they have identified a novel herpesvirus.

The need for further biological characterization to clarify the relationship between a novel pathogen and a disease is a wellrecognized problem in viruses studied primarily from molecular data (reviewed in Lipkin, 2010). For example, Bodewes et al. (2014) explicitly comment on the problem of establishing etiology in their report of a novel herpesvirus in harbor and gray seals, initially discovered in juvenile seals with ulcerative gingivitis but also found to be present in a large proportion of healthy controls.

### **Methodology Influences Success**

Metagenomic and deep-sequence analysis of a wide range of eukaryotic organisms over the last decade has shown that herpesviruses are more diverse and ubiquitous than previously imagined. Successful identification of novel herpesviruses or related viruses will rely on the availability of high-quality DNA and RNA from suitable primary samples.

While samples collected non-invasively, such as feces or urine, are immediately appealing for virus discovery because they are easy to collect, they may not be the ideal source material for identifying novel herpesviruses. Metagenomic analysis of bats found most herpesvirus-related sequence reads in throat swab or saliva samples. The primate sequence data in which novel herpesviruses were identified in public sequence repositories were initially sampled from blood and liver. There is an obvious correlation with the cellular tropism of individual herpesvirus species, their abundance in a particular tissue, and the likelihood of discovering a herpesvirus. Herpesviruses with a tropism for lymphocytes may be most readily detectable in blood samples or lymph node biopsies, for example. The physiological and environmental stresses affecting the host may also be an important variable in novel herpesvirus detection.

Metagenomics is an important tool in characterizing herpesviruses previously only known of from sequencing of short PCR amplicons. The success of metagenomics methods improves as the proportion of viral reads within a sample increases. Enriching for total viral material within a sample can be achieved by physical methods early in the sample preparation process, physically reducing the amount of contaminating non-viral nucleic acid using combinations of centrifugation, filtration and nuclease treatments (reviewed in Hall et al., 2014). Further enrichment is then possible, either sequence-independent enrichment using random hexamers or targeted enrichment (such as Sure Select) of specific viruses. All of these approaches may be needed to assemble whole genomes of herpesviruses from uncultured primary samples (Depledge et al., 2011). Targeted enrichment of herpesviruses relies on whole genome sequences being available for the virus of interest. Deep sequencing data from metagenomic studies can provide the sequences needed to design an initial set of PCR primers, as with phocine herpesvirus 7 (Bodewes et al., 2014), or target enrichment baits so that further sequence information can be collected. Metagenomics has also identified novel viruses from many other families, including coronaviruses (Schurch et al., 2014), papillomaviruses (Canuti et al., 2014), and tornoviruses (Ng et al., 2009), to give just a few examples.

Our closest living relatives, chimpanzees, bonobos and gorillas, are a rich source of human herpesvirus-like agents; whole-genome sequences of these viruses would tell us about our shared evolution with these ancient and successful pathogens. The discovery of two previously unknown primate herpesviruses in public sequence repositories (Aswad and Katzourakis, 2014) makes clear the possibilities for identifying further novel primate herpesviruses through metagenomics in the near future.

There is also the possibility of discovering further human herpesviruses. Humans are unusual among great apes because they have only one rhadinovirus species (KSHV), not two. The presence of three highly divergent forms of KSHV gene K15, including a form found only in southern Africa, has been suggested by some authors (Hayward and Zong, 2007) to indicate that there may be (or may have been) an unknown second KSHV-like rhadinovirus in humans. Deep sequencing and metagenomic analysis of primary tissue samples has been successful in identifying other primate herpesviruses, and as our knowledge of herpesvirus sequence diversity increases, additional human herpesviruses may be identified.

### **References**

Metagenomics is a promising technique that provides a further tool to identify novel herpesviruses, reducing our reliance on PCR for herpesvirus discovery while also improving our ability to design PCR primers or baits for sequencing studies. Furthermore, metagenomics is increasing our knowledge of the diversity and evolutionary history of herpesviruses; and emphasizing the importance of studying primary samples rather than cell lines or cultured virus isolates if we wish to discover further novel herpesviruses.

### **Acknowledgments**

CJH was funded by GOSH/UCL BRC. JB receives funding from NIHR (www.nihr.co.uk) and UCLH/UCL BRC (http://www. uclhospitals.brc.nihr.ac.uk/).

Alaska, through iterative assembly of taxonomic units. *PLoS ONE* 9:e105227. doi: 10.1371/journal.pone.0105227

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Copyright © 2015 Houldcroft and Breuer. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

## High temporal and spatial diversity in marine RNA viruses implies that they have an important role in mortality and structuring plankton communities

### *Julia A. Gustavsen1, Danielle M. Winget 1†, Xi Tian2 and Curtis A. Suttle1,3,4\**

*1 Department of Earth, Ocean and Atmospheric Sciences, University of British Columbia, Vancouver, BC, Canada*

*2 Bioinformatics Graduate Program, Faculty of Science, University of British Columbia, Vancouver, BC, Canada*

*3 Departments of Botany, and Microbiology & Immunology, University of British Columbia, Vancouver, BC, Canada*

*4 Canadian Institute for Advanced Research, Toronto, ON, Canada*

#### *Edited by:*

*Bas E. Dutilh, Radboud University Medical Center, Netherlands*

#### *Reviewed by:*

*Alexander Culley, Université Laval, Canada Jennifer R. Brum, University of Arizona, USA*

#### *\*Correspondence:*

*Curtis A. Suttle, Department of Earth, Ocean and Atmospheric Sciences, University of British Columbia, 2020 - 2207 Main Mall, Vancouver, BC V6T 1Z4, Canada e-mail<EMAIL_ADDRESS>
#### *†Present address:*

*Danielle M. Winget, Department of Biology, Seattle Pacific University, Seattle, WA, USA*

Viruses in the order *Picornavirales* infect eukaryotes, and are widely distributed in coastal waters. Amplicon deep-sequencing of the RNA dependent RNA polymerase (RdRp) revealed diverse and highly uneven communities of picorna-like viruses in the coastal waters of British Columbia (BC), Canada. Almost 300 000 pyrosequence reads revealed 145 operational taxonomic units (OTUs) based on 95% sequence similarity at the amino-acid level. Each sample had between 24 and 71 OTUs and there was little overlap among samples. Phylogenetic analysis revealed that some clades of OTUs were only found at one site; whereas, other clades included OTUs from all sites. Since most of these OTUs are likely from viruses that infect eukaryotic phytoplankton, and viral isolates infecting phytoplankton are strain-specific; each OTU probably arose from the lysis of a specific phytoplankton taxon. Moreover, the patchiness in OTU distribution, and the high turnover of viruses in the mixed layer, implies continuous infection and lysis by RNA viruses of a diverse array of eukaryotic phytoplankton taxa. Hence, these viruses are likely important elements structuring the phytoplankton community, and play a significant role in nutrient cycling and energy transfer.

**Keywords: RNA viruses, phytoplankton mortality, viral ecology, pyrosequencing, viral diversity,** *Picornavirales*

### **INTRODUCTION**

Viruses are highly abundant and widespread in the oceans (Bergh et al., 1989; Suttle, 2005). Beyond their impacts on host mortality, viruses are significant mediators of biogeochemical processes, horizontal gene transfer, and host community diversity in the oceans (Fuhrman, 1999; Wilhelm and Suttle, 1999; Suttle, 2005). Marine viruses are important pathogens of phytoplankton (Brussaard, 2004), and have been implicated in the termination of blooms (Nagasaki et al., 1994; Schroeder et al., 2003) and with succession in phytoplankton communities (Mühling et al., 2005). Viruses have been characterized that infect a wide variety of phytoplankton such as haptophytes (Bratbak et al., 1993), prasinophytes (Mayer and Taylor, 1979; Brussaard, 2004; Brussaard et al., 2004; Derelle et al., 2008), chlorophytes (Van Etten et al., 1981), diatoms (Shirai et al., 2008), and dinoflagellates (Tomaru et al., 2004). Viruses infecting eukaryotic phytoplankton generally have very narrow host ranges (Short, 2012).

Viruses infecting marine phytoplankton have genomes comprised of double-stranded DNA (dsDNA), single-stranded DNA (ssDNA), double-stranded RNA (dsRNA), and single-stranded RNA (ssRNA) (as reviewed in Short, 2012). Their genomes and particle sizes range from very large dsDNA viruses in the *Phycodnaviridae* to very small ssDNA and ssRNA viruses belonging to the genus *Bacilladnavirus* and order *Picornavirales*, respectively. The order *Picornavirales* is comprised of positivesense, ssRNA viruses that infect eukaryotes (Le Gall et al., 2008), including ecologically important marine protists. These viruses are small (25–35 nm), icosahedral, and have a conserved genomic organization that includes a replication area comprised of a type III helicase, a 3C-like proteinase, and a type I RNA dependent RNA polymerase (Sanfaçon et al., 2009). Isolates in the *Picornavirales* that are pathogens of marine protists infect a wide diversity of hosts including the bloom-forming raphidophyte *Heterosigma akashiwo* (Tai et al., 2003) (viral family *Marnaviridae*), the thraustochytrid *Aurantiochytrium s*p. (Takao et al., 2006; Yokoyama and Honda, 2007) (proposed viral genus *Labyrnavirus*) and the cosmopolitan diatoms *Rhizosolenia setigera* (Nagasaki et al., 2004) and *Chaetoceros socialis* (Tomaru et al., 2009) (proposed viral genus *Bacillarnavirus*). Viruses in the *Picornavirales* appear to be common and widely distributed in coastal waters (Culley et al., 2003; Culley and Steward, 2007).

Metagenomic and targeted gene studies are uncovering the diversity of marine RNA viruses. For example, phylogenetic analysis of RNA-dependent RNA polymerase (RdRp) sequences from seawater samples supports a monophyletic marine group within the *Picornavirales* (Culley et al., 2003, 2014; Culley and Steward, 2007; Tomaru et al., 2009) and several divergent clades within this marine group (Culley et al., 2003, 2014; Culley and Steward, 2007). Additionally, metagenomic analyses reveal that there are numerous sequences from aquatic RNA viruses that cannot be assigned to known taxa (Culley et al., 2006, 2014; Djikeng et al., 2009; Steward et al., 2012). Despite the high diversity of marine RNA viruses (Lang et al., 2009), the spatial and temporal distribution of different phylogenetic groups remains unreported, although there is evidence that the taxonomic structure of marine RNA viral communities is highly uneven. For example, in one sample from a metagenomic study from the coastal waters of British Columbia, 59% of the reads assembled into a single contig, while in a second sample 66% of the reads fell into four contigs, with most falling into two genotypes (Culley et al., 2006). However, with only a few 100 reads in total from the two samples, the coverage of the communities was low. Similarly, RNA viral metagenomic data from a freshwater lake (Djikeng et al., 2009) showed little identical sequence overlap among communities, although there was broad taxonomic similarity over time within a location.

Ecological questions about the distribution of marine viruses over time and space have been examined more extensively in bacteriophages, particularly those infecting cyanobacteria. For example, some data reveal no clear patterns of biogeography in cyanophage isolates locally (Clasen et al., 2013), regionally (Jameson et al., 2011) or more globally (Huang et al., 2010). Other data have shown patterns at a regional scale (Marston et al., 2013), although communities in basins that were connected were most similar and those that were separated by land or current boundaries were the least similar. Other data for marine bacteriophages have shown temporal variability (Chen et al., 2009; Wang et al., 2011; Chow and Fuhrman, 2012; Clasen et al., 2013; Marston et al., 2013). If the dynamics of marine bacteriophages and marine RNA viruses are similar, some RNA viral taxa will persist temporally and spatially, while other taxa will be detected sporadically. To test this hypothesis, we examined two samples, taken 5 months apart at the same location, and three samples taken within hours of each other, but 20 km apart in the same coastal basin.

We used high-throughput 454 pyrosequencing to obtain deep coverage of RdRp amplicon sequences and compare the richness of viruses in the *Picornavirales* among samples from the coastal waters of British Columbia, Canada. The results revealed a phylogenetically diverse and spatially variable community of viruses, suggesting that taxon-specific lytic events are important in shaping the phytoplankton community.

### **MATERIALS AND METHODS SAMPLING LOCATIONS**

To assess viral communities from different coastal habitats, we collected samples from three sites in the Strait of Georgia (49◦14.926N 123◦35.682W, 49◦17.890N 123◦43.650W, and 49◦23.890N 123◦59.706W), and from Jericho Pier (49◦16 36.73N, 123◦12 05.41W) in British Columbia, Canada) (**Figure 1**). The Strait of Georgia (SOG) is an estuarine-influenced basin that is on average 22 km across, 222 km long and 150 m deep. The upper 50 m of SOG is where most of the variability in physical and

British Columbia, Canada. Jericho Pier was sampled in summer (JP-S) and fall (JP-F).

chemical parameters occurs. Jericho Pier (JP) is adjacent to the shoreline, in a well-mixed location with mixed semidiurnal tides.

#### **SAMPLE COLLECTION**

On 28 July 2010, a rosette equipped with a Seabird SBE 25 CTD (equipped with Seabird SBE 43 dissolved oxygen sensor, WET Labs WETStar fluorometer, WET Labs C-Star transmissometer, Biospherical Instruments QSP2200PD PAR sensor) and General Oceanics GO-FLO bottles was used to collect 12 L of water from five depths between 2 and 16 m at each of three stations in SOG. The five depths from each station were combined to obtain three integrated water samples. The depths were selected to attempt to encompass the viral diversity at each station, and consisted of a near-surface sample, a sample from the isothermal zone below the mixed layer, and three samples spanning the chlorophyll maximum.

At JP, 60 L of water was pumped from the 1-m depth on 10 July and 12 October 2010. Salinity and temperature at JP were measured using a YSI probe (Yellow Springs, Ohio, USA). For all samples, the water was filtered through 142-mm diameter, 1.2-µm nominal pore-size glass-fiber (GC50 Advantec MFS, Dublin, CA., USA) and 0.22-µm pore-size polyvinyldine (Millipore Bedford, MA, USA) filters. The viral size fraction in the filtrate was concentrated to ∼500 mL (viral concentrate) using tangential flow ultrafiltration using a 30 kDa MW prep-scale Spiral Wound TFF-6 cartridge (Millipore) (Suttle et al., 1991).

Phosphate, silicate and nitrate+nitrite concentrations were determined in duplicate 15-mL seawater samples filtered through 0.45-µm pore-size HA filters (Millipore) and stored at −20◦C until air-segmented continuous-flow analysis on a AutoAnalyzer 3 (Bran & Luebbe, Norderstedt Germany). Chlorophyll *a* (Chl *a*) was determined in triplicate by filtering 100 mL of seawater onto 0.45µm pore-size HA filters (Millipore), and storing the filters in the dark at −20◦C until acetone extraction and then analyzed fluorometrically (Parsons et al., 1984). The average and standard error of the replicates was calculated for each sample.

### **NUCLEIC-ACID EXTRACTION AND PCR**

The viral concentrate was filtered twice through 0.22-µm pore-size Durapore PVDF filters (Millipore) in a sterile Sterivex filter unit (Millipore). The filtrate, containing virus-sized particles, was pelleted by ultracentrifugation (Beckman-Coulter, Brea, California, USA) in a SW40 rotor at 108 000 *g* for 5 h at 12◦C. The pellet was resuspended overnight in 100µL of supernatant at 4◦C. To digest free DNA, the pellets were incubated with 1U/µl DNAse with a final concentration 5 mM MgCl2 for 3 h at room temperature. Nucleic acids were extracted using a Qiamp Viral Minelute spin kit (Qiagen, Hilden, Germany) according to the manufacturer's directions. To remove DNA, the extracted viral pellets were digested with DNase 1 (amplification grade) (Invitrogen, Carlsbad, California, USA) and the reaction was terminated by adding 2.5 mM EDTA (final concentration) and incubating for 10 min at 65◦C. Complementary DNA (cDNA) was generated using Superscript III reverse transcriptase (Invitrogen) with random hexamers (50 ng/µl) as per the manufacturer. PCR was performed using primer set MPL-2 for a targeted set of the marine picornavirus-like RdRp (Culley and Steward, 2007). Each reaction mixture (final volume, 50µl) consisted of 50 ng of cDNA, 1× (final concentration) PCR buffer (Invitrogen), 2 mM MgCl2, 0.2 mM of each deoxynucleoside triphosphate (Bioline, London, UK), 1µM of each primer, and 1 U Platinum Taq DNA polymerase. The reaction was run in a PCR Express thermocycler (Hybaid, Ashford, UK) with the following conditions: 94◦C for 75 s, followed by 40 cycles of denaturation at 94◦C for 45 s, annealing at 43◦C for 45 s, and extension at 72◦C for 60 s and a final extension step of 9 min at 72◦C. PCR products were cleaned using the Minelute PCR purification kit (Qiagen).

### **LIBRARY PREP AND PYROSEQUENCING**

Libraries for each site were prepared for sequencing using NEBNext DNA Library Prep for 454 kit (New England Biolabs, Ipswich, Massachusetts, USA) following the manufacturer's directions, and using Ampure beads (Beckman-Coulter) for size selection and purification using a bead ratio of 0.8:1 beads:library. PCR amplicons were barcoded and sent for 454 Titanium pyrosequencing (Roche, Basel, Switzerland) at the Génome Québec Innovation Centre at the McGill University (Montreal, QC, Canada).

#### **SEQUENCE ANALYSIS**

Sequencing reads were quality trimmed using length settings between 100 and 600 bp, with a maximum of 3 primer mismatches for the specific primer, and denoised using the denoiser algorithm in QIIME (version 1.7) with default settings for Titanium data (Caporaso et al., 2010; Reeder and Knight, 2010). Sequences were checked for chimeras using UCHIME (Edgar, 2010) against a nucleotide database of RdRp sequences built using NCBI-BLAST+ by retrieving nucleotide sequences of the RdRp from *Picornavirales* viral isolates (accessed: 8 August 2012), and also using the *denovo* chimera check in UCHIME. The overlap of results from these two methods was defined as chimeric sequences, although none were found in this study. Non-chimeric sequences were queried using BLASTx (Altschul et al., 1990) with an e-value of 1e−3 against the database of RdRp viral isolates. All sequences with hits were retained and all sequences with no hits were then queried against the non-redundant (nr) Genbank database (Benson et al., 2007) using BLASTx with an e-value of 1e−5. All the sequences identified as contaminants or as unknown with only 1 read were removed. Remaining sequences were translated to amino acids using FragGeneScan with the 454\_10 training option (Rho et al., 2010). Sequences were grouped into operational taxonomic units (OTUs) using UCLUST at a range of similarities from 50 to 100% (Figure S1) using the original seed sequences (centroids) as the output (Edgar, 2010). A similarity of 95% was chosen for this analysis for the following two main reasons: (1) When the NCBI conserved domain alignment for the RdRp region (all *Picornavirales)* was analyzed for percent similarity, the only sequences that displayed greater than 95% similarity in this region were strains of the same virus. (2) The sequences from the control libraries (consisting of 1 clone) clustered into 1 sequence at this percentage (see Supplemental methods). Thus clustering at 95% similarity was a way to use biological and sequence-based information to inform our choice of cut-off to collapse strain level variation and as a conservative approach to avoid variation that may be present because of the sequencing platform. Control sequences obtained by cloning and Sanger sequencing (see Supplemental Methods) were used to verify the sequence processing methodology. Raw and processed sequence data were deposited in the NCBI BioProject database ID: PRJNA267690.

#### **PHYLOGENETIC ANALYSIS**

All OTUs with less than 5 reads were removed and the remaining OTUs were aligned using profile alignment in Muscle (Edgar, 2004) to seed alignments of viral RdRps from the NCBI Conserved Domain Database (Marchler-Bauer et al., 2011). Sequences from other environmental surveys were clustered in the same manner as the reads in this study (using UCLUST at 95%).  Alignments were masked using trimAl with the automatic heuristic (Capella-Gutierrez et al., 2009) and edited manually. ProtTest 3.2 was used for amino-acid model selection (Darriba et al., 2011) before building the initial maximum likelihood phylogenetic tree using FastTree (Price et al., 2010). Final trees were done with RA×ML using sequences belonging to viruses in the *Sequiviridae* as the outgroup, and the BLOSUM62 amino-acid model with 100 bootstraps (Stamatakis et al., 2008). The tree was visualized in R (R Core Team, 2014) using the ape package (Paradis, 2012) and edited in Figtree (http://tree*.*bio*.*ed*.*ac*.*uk/software/figtree/).

#### **STATISTICAL ANALYSIS**

Generation of rarefaction curves by random resampling of OTU abundances was performed using the vegan package (Oksanen et al., 2013) in R (R Core Team, 2014). Relative abundances were normalized by randomly resampling 10 000 times using vegan (Oksanen et al., 2013), normalizing to the library with the lowest number of reads and then taking the median. Rank-abundance curves were generated with ggplot2 (Wickham, 2009) using the normalized OTU relative abundances. Scripts used in this project are available as part of QIIME and custom user scripts used to process the data are available here: http://github.com/jooolia/RdRp\_454\_amplicons\_Jericho\_ and\_SOG; doi: 10.5281/zenodo.12509

### **RESULTS**

#### **ENVIRONMENTAL PARAMETERS**

The environmental parameters ranged widely among samples (**Figure 2**). Chlorophyll *a* values were lowest at Jericho Pier Fall (JP-F) at 0.16µg L−1 (± 0.04) and highest at the Strait of Georgia Station 2 (SOG-2) at 3.2µg L−1 (± 0.4). Silicate values for all samples were similar (range of 25.5 to 37.8 µM), except for Jericho Pier Summer (JP-S) when silicate was lower at 6.2µM (± 0.01). Phosphate ranged between 0.90 and 1.3µM at JP-F, SOG-1, and SOG-4, but was lower at SOG-2 (0.58µM (± 0.03)) and lowest at JP-S (0.06µM (± 0.04)). Nitrate + nitrite values were more variable than the other nutrients and ranged from 1.35µM (± 0.75) at JP-S to 14.6µM (± 0.04) at JP-F. The SOG sites were highly stratified with SOG-4 being the most stratified with a calculated mixed-layer depth of 2 m, while SOG-1 and SOG-2 were similar with a mixed-layer depth of 6 m (**Table 1**).

#### **ANALYSIS OF RdRp SEQUENCES**

After quality filtering to remove homopolymers and contaminating reads, 300 180 reads were recovered from the 5 libraries of RdRp amplicons. At all sites the rarefaction curves plateaued indicating that the depth of sampling was adequate to assess the communities (**Figure 3**). From these reads, 265 unique OTUs (at 95% similarity) were identified, including 108 singletons. For further analysis OTUs were excluded that did not contain recognizable RdRp motifs (Koonin, 1991; Le Gall et al., 2008), generally did not align well with other RdRp sequences and those that were not present in any sample after normalization.

Using the above criteria there were 145 OTUs identified in all samples, and between 24 to 71 OTUs per site. The Jericho Pier samples had 116 OTUs of which only 10 (8.6%) were shared between sampling times (**Figure 4**). JP-S had the highest richness with 71 OTUs, of which 59 (83 %) were unique, while JP-F had the second highest richness (49 OTUs), of which 45 (92 %) were unique. The SOG sites together had 64 OTUs, none of which were shared among all sites. SOG-1 and SOG-2 had the lowest number of OTUs (24). SOG-1 had only three OTUs which were unique. However, 21 (33%) were shared between SOG-1 and SOG-4, and 6 (9%) between SOG-2 and SOG-4. The majority of OTUs (75%) from SOG-2 were unique, whereas most OTUs from the other SOG sites were shared with other sites (87% for SOG-1 and 63% for SOG-4).

Rank abundance curves of the viral OTUs showed that at each site most sequences were assigned to only a few OTUs (**Figure 5**).

**FIGURE 2 | Environmental parameters. (A)** Chlorophyll *a* with standard error of the mean from triplicates. **(B)** Silicate with standard error of the mean from duplicates. **(C)** Phosphate with standard error of the mean from duplicates. **(D)** Nitrate+ nitrite with standard error of the mean from

duplicates. **(E)** Temperature with each point as one GO-FLO bottle (SOG samples) or from total seawater sample (Jericho samples). **(F)** Salinity with each point representing a seawater sample as one GO-FLO bottle (SOG samples) or from total seawater sample (Jericho samples).

*\*after quality filtering and matching to the RdRp primer set.*

JP-S had the highest richness but the shallowest slope of these curves, demonstrating more evenness in the abundance of OTUs than at the other sites. SOG-4 and JP-S had similar rank abundance curves that were much shallower that those of SOG-1 and SOG-2 (**Figure 5**).

The OTUs that were observed in more than 5 reads were placed in phylogenetic context using a maximum likelihood RA×ML tree (Stamatakis et al., 2008) with sequences from previous RdRp gene surveys and isolated viruses (**Figure 6**). OTUs from this study fell within a well-supported clade that includes all the marine isolates belonging to the *Picornavirales*. Within this group there was a well-supported divide between OTUs grouping in the *Marnaviridae* clade and those grouping with sequences from viruses infecting diatoms and a thraustochytrid. The overall tree topology is not well supported, although there are a number of well-supported clades containing OTUs from this study and other environmental sequences.

The *Marnaviridae* clade had the greatest number of OTUs (10) associated with it; whereas, very few OTUs (only OTUs 89, 107, 75, and 120) from this study were assigned to clades primarily from Hawaii (Culley and Steward, 2007; Culley et al., 2014). No clade contained OTUs from all sites. The Jericho Pier samples were the most phylogenetically diverse (**Figure 6**, **Table 2**), and contained OTUs (e.g., OTUs 6, 7, 35, 31, 84, 47, 14, 4, 34, 39, 44, 23, 8, 20) that fell into clades that did not contain OTUs from any of the SOG samples. Some clades contained OTUs from both JP-S and JP-F samples; however, many OTUs within the clades were unique to one Jericho Pier sample. Phylogenetic diversity differed among samples, except for OTUs from the two SOG sites with deeper mixed layers, some of which were present in different clades resulting in similar phylogenetic diversity (**Table 2**).

The Strait of Georgia (SOG) sites were sampled within hours of each other, and the water at each site was pooled from multiple depths above, below, and across the chlorophyll maximum. One of the most striking differences among sites was that SOG-1 and SOG -2 had mixed layer depths of 6 m; whereas SOG-4 had a mixed-layer depth of 2 m, and much higher richness and phylogenetic diversity. Sites SOG-1 and SOG-2 had the lowest phylogenetic diversity (**Table 2**). All the OTUs found at SOG-1 (33, 12, 16, 9, 10, 27, 29) were within the *Marnaviridae* clade; similarly, all OTUs (5, 82, 2, 1) from SOG-2 were within one distantly related clade. In both cases OTUs from these clades occurred at SOG-4. SOG-2 did not have the high numbers of HaRNAV-related viruses that were found in all other samples.

#### **DISCUSSION**

Pyrosequencing of RdRp gene fragments from coastal samples uncovered much greater genetic diversity than in previous gene surveys (Culley et al., 2003, 2014; Culley and Steward, 2007) and revealed many previously unknown taxonomic groups within the *Picornavirales*. As well, striking differences in the taxonomic richness among samples implies that these viruses infect a wide variety of eukaryotic plankton, but that the mortality imposed on some taxa is highly variable across space and time. Other taxonomic groups within the *Picornavirales* were more widespread, suggesting that infection of some planktonic taxa is more widespread and persistent. These results and their implications are discussed in detail below.

#### **EXPANDING THE KNOWN DIVERSITY OF** *PICORNAVIRALES*

The high depth of sequencing and limited diversity in each library (**Figure 3**) gives high confidence that the population structure of RdRp amplicons in each sample has been well characterized (Kemp and Aller, 2004). Although some sequences were closely related to those found in previous studies (Culley et al., 2003; Culley and Steward, 2007) (**Figure 6**), many OTUs formed new clades. Many OTUs were related to *Heterosigma akashiwo* RNA virus (HaRNAV) that infects the toxic bloom-forming raphidophyte *Heterosigma akashiwo* (Tai et al., 2003). HaRNAV is the

type virus of the family *Marnaviridae* (Lang et al., 2004); it has a genome of about 9.1 kb and a high burst size as indicated by the large crystalline arrays of particles in the cytoplasm of infected cells (Tai et al., 2003). HaRNAV was isolated from coastal waters in British Columbia, Canada (Tai et al., 2003) and can remain infectious for many years in sediments (Lawrence and Suttle, 2004). Interestingly, HaRNAV was isolated from the same area as the present study, and appeared ancestral to many of the recovered sequences based on the phylogeny. For example, OTU 0 was most abundant (18,034 reads after rarefaction) and clustered closely with HaRNAV, although the sequence was only 76.4% similar at the amino-acid level. However, other OTUs in the cluster ranged between 55 and 79% similar to HaRNAV, which is low compared to amino-acid similarities of other RNA viruses within a family that usually have greater than 90% aa similarity (Ng et al., 2012).

#### **DISTINCT COMMUNITIES OCCURRED IN DIFFERENT SEASONS AT THE SAME LOCATION**

While only 8.6% of the OTUs from Jericho Pier were shared between dates, both samples had similar evenness, although the summer sample had greater richness (**Figure 5**, **Table 2**). The small overlap in OTUs between sampling dates is not surprising given the very different conditions between July and October (**Figure 1**), and the dynamic nature of planktonic communities in response to environmental changes. At the same location, dsDNA viruses belonging to the *Phycodnaviridae*, which infect eukaryotic phytoplankton, varied seasonally based on fingerprint analyses of DNA polymerase gene fragments using denaturing gradient gel electrophoresis; however, some OTUs persisted for extended periods (Short and Suttle, 2003). Similarly, the composition of other aquatic viral communities has been shown to be dynamic although some OTUs persist (Djikeng et al., 2009; Rodriguez-Brito et al., 2010), and in some cases have repeatable seasonal patterns (Chow and Fuhrman, 2012; Clasen et al., 2013; Marston et al., 2013). With only single samples from summer and fall, inferences about dynamics cannot be made from our data.

One of the few taxonomic groups that occurred in both the summer (JP-S) and fall (JP-F) samples from Jericho Pier was related to HaRNAV (**Figure 6**). There was greater diversity of OTUs in this clade in JP-F, even though JP-S had higher richness and higher phylogenetic diversity overall. This is unlike bacterial and phytoplankton communities that tend to be more diverse in winter (Zingone et al., 2009; Ladau et al., 2013). However, the

RdRp primers target a specific subset of the viral community that might not reflect the overall taxonomic diversity.

*Picornavirales* sequences. Bootstrap values above 65% are labeled. Adjacent

Based on genome organization and sequence identity, RNA viruses that infect diatoms have been assigned to the genus *Bacillarnavirus*, that includes *Rhizosolenia setigera* RNA virus (RsRNAV) (Nagasaki et al., 2004)*, Chaetoceros tenuissimus* Meunier RNA virus (CtenRNAV) (Shirai et al., 2008) and *Chaetoceros socialis* f. radians RNA virus (CsfrRNAV) (Tomaru et al., 2009). In the JP-F sample, the most relatively abundant cluster grouped with RsRNAV that infects the marine diatom *Rhizosolenia setigera* (Nagasaki et al., 2004). This corresponded with the highest levels of nitrate + nitrite, which is often associated with high diatom abundances (Zingone et al., 2009); hence, these OTUs are likely associated with viruses infecting diatoms.

#### **DISTINCT COMMUNITIES OCCURRED AT GEOGRAPHICALLY PROXIMATE SITES**

the sample with the lowest number of reads.

Areas of higher habitat diversity, such as stratified water layers, generally have higher biological richness (Klopfer and MacArthur, 1960; Chesson, 2000), and this is consistent with

**Table 2 | Phylogenetic diversity, species richness.**

*Phylogenetic diversity (PD) is calculated as in Faith (1992). OTUs were the same as used in the construction of the phylogenetic tree and must have included 5 or more reads.*

the much higher richness and phylogenetic diversity found at SOG-4, which was the most stratified site and included the most abundant OTUs from SOG-1 and SOG-2. Most OTUs from SOG-1 clustered in the *Marnaviridae* clade, while most SOG-2 OTUs clustered in a phylogenetically distant clade. Given that we have used very conservative clustering, and that dsDNA viruses infecting phytoplankton are strain specific and have phylogenies that are congruent with their hosts (Clasen and Suttle, 2009; Bellec et al., 2014), and that RNA viruses infecting diatoms and the dinoflagellate, *Heterocapsa circularisquama* (Nagasaki et al., 2005) are host-specific, it implies that closely related OTUs infect closely related taxa of phytoplankton. Hence, it suggests that the most abundant viruses at these three locations infect different species.

There are few clear patterns in the spatial distribution of viruses in marine waters where geographically distant sites are connected by currents and mixing. The best examples are for cyanophages. For instance, when looking at local variation in cyanophages isolated at sites in Southern New England, 72% of the viral OTUs were shared between at least 2 sites (Marston et al., 2013); however, between Bermuda and Southern New England only 2 OTUs overlapped and they comprised only 0.6% of the isolates. Yet, clear patterns of cyanophage OTU distribution by depth occurred in areas adjacent to the SOG when assessed using community fingerprinting (Frederickson et al., 2003). The biggest differences with depth occurred in stratified water in which some OTUs were present at all depths, while others were only present at specific depths, even though the samples were collected only meters apart (Frederickson et al., 2003). These viruses infect cyanobacteria, as opposed to the picorna-like viruses, which likely infect protistan plankton. Nonetheless, the factors governing the distribution of cyanobacterial and protistan hosts are likely similar; hence, different OTUs would be expected to occur in environments with different vertical structure (stratification) of the water column.

Rank abundance curves showed that SOG-1 and SOG-2 were the least even communities (**Figure 5**). Overall, at most sites four to five viral OTUs were most abundant (**Figure 5**) similar to other reports for aquatic viral communities in which a few viruses dominate, but most of the diversity comes from rarer viruses (Angly et al., 2006; Suttle, 2007). Our targeted approach showed that the picornavirus-like virus communities at SOG and JP were dominated by only a few genotypes, supporting previous metagenomic results showing that the OTU distributions of RNA viruses in SOG and JP were highly uneven with little overlap between sites (Culley et al., 2006).

#### **EACH OTU LIKELY REPRESENTS A SINGLE LYTIC EVENT**

Given that the hosts of marine *Picornavirales* isolates are protists, and that protists are the most abundant eukaryotes in the sea, it is likely that the majority of OTUs recovered in this study are from viruses that infect these unicellular marine eukaryotes. These eukaryotic communities are highly dynamic and change throughout the year based on environmental and biological factors (Larsen et al., 2004). Since viral infection is usually host specific, the diversity in marine viral communities is a reflection of the underlying diversity of the marine eukaryotic hosts. Moreover, viral propagation is dependent on host encounter rates and is proportional to host-cell abundance (Murray and Jackson, 1992); hence the most abundant taxa will be most likely to encounter and propagate a viral infection, giving the opportunity for rarer species to increase in abundance and promoting diversity (Thingstad, 2000; Winter et al., 2010). Since our study was not over time it is difficult to evaluate whether these data support the Bank model (Breitbart and Rohwer, 2005), however, some taxa were found at one site, but not a similar nearby site, thus these taxa could be present at background levels at some sites and more abundant in others.

It is probable that the most abundant OTUs in these data are from recent lysis of host taxa. An error rate for replication of RNA viruses of about 1 bp mutation per generation (9000 bp genome × 0.0001 error rate per base pair = 1 bp; Holmes, 2009), and a lower-end burst-size estimate of 1000 particles for marine viruses in the *Picornavirales* that infect protists (Lang et al., 2009), would produce about 1000 different genomes from each lysed cell. For the amplified 500 bp RdRp gene fragment there is a 0.00056% chance of an error in 1 generation, assuming that mutations are distributed evenly in the genome (Sanjuan et al., 2010; Combe and Sanjuán, 2014). Consequently, even with the relatively high error rates of RNA replication, when grouped at 95% similarity at the amino acid level, all of the sequences from a lytic event should fall within a single OTU. The halflife for decay of viral infectivity and particles in the surface mixed layer is typically a few hours (Heldal and Bratbak, 1991; Suttle and Chen, 1992; Noble and Fuhrman, 1997; Bettarel et al., 2009); thus the recovered viral OTUs were likely from recent lytic events. Furthermore, considering the specificity of viruses infecting protists (Short, 2012), each OTU probably stems from viruses infecting a single host taxon. Thus, these data imply that infection of marine protists by viruses in the *Picornavirales* is not only pervasive, but likely involves a wide diversity of host taxa; hence, these viruses are likely important structuring elements for phytoplankton communities that influence nutrient cycling and energy flow.

#### **AMPLICON DEEP SEQUENCING AS AN APPROACH FOR ESTIMATING VIRAL DIVERSITY**

Amplicon deep sequencing is a sensitive and high-resolution approach for examining microbial community dynamics over time and space (Caporaso et al., 2011; Gobet et al., 2012; Gibbons et al., 2013). Careful quality trimming of sequences and removal of singletons is essential for reliable results (Zhou et al., 2011) since errors in sequences will inflate estimates of diversity. With careful data processing and analysis, amplicon deep sequencing is as accurate for assessing community composition and diversity as cloning and Sanger sequencing (Amend et al., 2010), but with much greater depth of coverage of the community.

There are potential biases associated with reverse transcription with random hexamers (which can decrease yield and could inflate diversity) (Zhang and Byrne, 1999), template amplification by PCR (Lee et al., 2012) and with using highly degenerate primers that target a specific part of the community containing many different templates (Culley and Steward, 2007). A danger of the high cycle number can be diversity overestimates which can come from the increasing number of chimeric sequences produced with greater cycle number (Qiu et al., 2001). The sequences were processed with caution considering the high number of PCR cycles employed in this study. Chimera checking *denovo* was used to look for chimeric sequences originating from two higher abundance reads, and reference-based chimera checking was used a database of RdRps from isolated viruses to correct for this potential error. In addition, a conservative cut-off was used of only OTUs comprising more than 5 reads that aligned to the conserved domain alignment.

Read abundance of OTUs can be considered semi-quantitative and good for comparisons of richness and diversity among samples (but not for absolute counts of genes) (Amend et al., 2010; Pinto and Raskin, 2012; Ibarbalz et al., 2014). Moreover, by using control sequences obtained by cloning and Sanger sequencing alongside pyrosequenced libraries containing the same sequence (Supplemental Methods, Figure S2) we verified that amplicon deep sequencing and our sequence processing methodology recovered accurate environmental viral sequences and non-inflated estimates of richness like in studies for bacterial amplicons (Sogin et al., 2006; Huse et al., 2008; Kirchman et al., 2010; Caporaso et al., 2011) and clinical viral studies (Romano et al., 2013; Watson et al., 2013).

### **CONCLUSION**

Amplicon deep sequencing of RdRp gene fragments using 454 pyrosequencing revealed the richness and population structure of marine *Picornavirales* in five coastal samples. The known diversity of viruses in this group was greatly increased with 145 OTUs that differed by at least 5% at the amino-acid level. There were between 24 and 71 OTUs in each sample, with distinct patterns of OTU distribution, richness and diversity among samples. There was little overlap between viral OTUs collected at the same site in summer and fall, and among samples collected 20 km apart on the same day. The high temporal and spatial diversity in RdRp sequences is consistent with viral communities that turnover rapidly, and episodic infection of a wide diversity of protistan hosts. The low overlap in OTUs and phylogenetic diversity among samples implies a dynamic landscape of viral infection and supports the idea that marine picorna-like viruses are important pathogens of marine protists that have an important role in structuring marine planktonic communities, and in nutrient cycling and energy transfer among trophic levels. Ultimately, further study is needed to disentangle the temporal and spatial drivers of these communities.

### **ACKNOWLEDGMENTS**

Thanks to R. Adelshin, A. M. Chan, C. Charlesworth, C. Chénard, R. Cruz, J. Finke, J. M. Labonté, J. Li, T. Nelson, J. P. Payet, E. J. Shelford, M. Vlok, R. A. White III, Q. Zhang, for help with sampling, and sample processing, the crew of the CGCS Vector for facilitating sample collection, and C. Payne for help with sampling, instrumentation and data acquisition, as well as processing of nutrient samples. M. Vlok and R. A. White III contributed to library preparation. Thanks to C. Chénard, C. Chow, J. Finke, T. Heger and E. J. Shelford for insightful comments on the manuscript. Funding was provided by NSERC Discovery and ship-time grants and from the Tula Foundation (Curtis A. Suttle), and NSERC PGS-M and PGS-D scholarships (Julia A. Gustavsen), and a UBC 4-Year Fellowship (Julia A. Gustavsen).

### **SUPPLEMENTARY MATERIAL**

The Supplementary Material for this article can be found online at: http://www.frontiersin.org/journal/10.3389/fmicb.2014.00703/ abstract

### **REFERENCES**

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 06 October 2014; accepted: 26 November 2014; published online: 15 December 2014.*

*Citation: Gustavsen JA, Winget DM, Tian X and Suttle CA (2014) High temporal and spatial diversity in marine RNA viruses implies that they have an important role in mortality and structuring plankton communities. Front. Microbiol. 5:703. doi: 10.3389/fmicb.2014.00703*

*This article was submitted to Virology, a section of the journal Frontiers in Microbiology.*

*Copyright © 2014 Gustavsen, Winget, Tian and Suttle. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# Previously unknown evolutionary groups dominate the ssDNA gokushoviruses in oxic and anoxic waters of a coastal marine environment

#### *Jessica M. Labonté1†, Steven J. Hallam1,2,3 and Curtis A. Suttle1,2,4,5\**

*1 Department of Microbiology and Immunology, University of British Columbia, Vancouver, BC, Canada, 2 Canadian Institute for Advanced Research, Toronto, ON, Canada, 3 Graduate Program in Bioinformatics, University of British Columbia, Vancouver, BC, Canada, 4 Department of Earth, Ocean and Atmospheric Sciences, University of British Columbia, Vancouver, BC, Canada, 5 Department of Botany, University of British Columbia, Vancouver, BC, Canada*

#### *Edited by:*

*Katrine L. Whiteson, University of California, Irvine, USA*

#### *Reviewed by:*

*Andrew David Millard, University of Warwick, UK Karyna Rosario, University of South Florida, USA*

#### *\*Correspondence:*

*Curtis A. Suttle, Department of Earth, Ocean and Atmospheric Sciences, University of British Columbia, 2178-2207 Main Mall, Vancouver, BC V6T 1Z4, Canada<EMAIL_ADDRESS>
#### *†Present address:*

*Jessica M. Labonté, Bigelow Laboratory for Ocean Sciences, East Boothbay, ME, USA*

#### *Specialty section:*

*This article was submitted to Virology, a section of the journal Frontiers in Microbiology*

> *Received: 27 October 2014 Accepted: 29 March 2015 Published: 22 April 2015*

#### *Citation:*

*Labonté JM, Hallam SJ and Suttle CA (2015) Previously unknown evolutionary groups dominate the ssDNA gokushoviruses in oxic and anoxic waters of a coastal marine environment. Front. Microbiol. 6:315. doi: 10.3389/fmicb.2015.00315* Metagenomic studies have revealed that ssDNA phages from the family *Microviridae* subfamily *Gokushovirinae* are widespread in aquatic ecosystems. It is hypothesized that gokushoviruses occupy specialized niches, resulting in differences among genotypes traversing water column gradients. Here, we use degenerate primers that amplify a fragment of the gene encoding the major capsid protein to examine the diversity of gokushoviruses in Saanich Inlet (SI), a seasonally anoxic fjord on the coast of Vancouver Island, BC, Canada. Amplicon sequencing of samples from the mixed oxic surface (10 m) and deeper anoxic (200 m) layers indicated a diverse assemblage of gokushoviruses, with greater richness at 10 m than 200 m. A comparison of amplicon sequences with sequences selected on the basis of RFLP patterns from eight surface samples collected over a 1-year period revealed that gokushovirus diversity was higher in spring and summer during stratification and lower in fall and winter after deep-water renewal, consistent with seasonal variability within gokushovirus populations. Our results provide persuasive evidence that, while specific gokushovirus genotypes may have a narrow host range, hosts for gokushoviruses in SI consist of a wide range of bacterial taxa. Indeed, phylogenetic analysis of clustered amplicons revealed at least five new phylogenetic groups of previously unknown sequences, with the most abundant group associated with viruses infecting SUP05, a ubiquitous and abundant member of marine oxygen minimum zones. Relatives of SUP05 dominate the anoxic SI waters where they drive coupled carbon, nitrogen, and sulfur transformations along the redoxline; thus, gokushoviruses are likely important mortality agents of these bacteria with concomittant influences on biogeochemical cycling in marine oxygen minimum zones.

#### Keywords: *Gokushovirinae*, host range, viral diversity, SUP05, oxygen minimum zones, PCR sequencing

## Introduction

Bacteriophages belonging to the *Microviridae* family consist of a ∼30 nm icosahedral capsid containing a positive-sense ssDNA molecule of 4.5–6.1 kb (King et al., 2012). Replication requires a minimum of two coat proteins (VP1 and VP2), a scaffolding protein (VP3), a replication protein (VP4), and a DNA packaging protein (VP5). Based on the phylogeny of the major capsid protein (MCP or VP1) from isolates, the *Microviridae* family is divided into two groups (Brentlinger et al., 2002); the *Microvirus* genus contains phages like phiX174 and G4 that infect *Escherichia coli* (Godson et al., 1978), while the *Gokushovirinae* subfamily includes those infecting parasitic bacteria such as *Chlamydia* [Chp1 (Storey et al., 1989), Chp2 (Liu et al., 2000; Everson et al., 2002), Chp3 (Garner et al., 2004)], *Bdellovibrio* [(phiMH2K; Brentlinger et al., 2002)], and *Spiroplasma* [(SpV4; Chipman et al., 1998)]. While it is commonly thought that *Microviridae* phages are strictly lytic (Liu et al., 2000; Garner et al., 2004; Salim et al., 2008), an *in silico* study found sequences with similar genome organization to gokushoviruses associated with *Bacteroidetes* from the human gut and mouth, which suggests that these phages can be temperate (Krupovic and Forterre, 2011). The temperate *Microviridae* phages are phylogenetically distinct from gokushoviruses, and have been assigned to the *Alpavirinae*, a proposed new sub-family within the *Microviridae* (Krupovic and Forterre, 2011). Recently, 81 microvirus genomes (including 42 gokushoviruses) were assembled from various environmental metegenomic data, identifying a new group, the *Pichovirinae*, which harbored a different genome organization of the conserved genes, indicating that microviruses display great diversity and may play an important role in many ecosystems (Roux et al., 2012b).

Marine gokushoviruses were first revealed in viral metagenomic data from the Strait of Georgia (SOG), Gulf of Mexico (GOM), and Sargasso Sea (SAR; Angly et al., 2006), and are among the most commonly recovered sequences from ssDNA phages in marine metagenomic data (Rosario and Breitbart, 2011). They were particularly abundant in the SAR, where 6% of the sequences were similar to the phage Chp1 that infects *Chlamydia psittaci* (Angly et al., 2006). This abundance of ssDNA sequences allowed for the assembly of two environmental *Gokushovirinae* genomes, with the help of PCR amplification (Tucker et al., 2011). A survey in the Atlantic Ocean showed a different depth distribution of these two genomes, consistent with alternative host-infection patterns. Sequences belonging to gokushoviruses have also been found in marine (Labonté and Suttle, 2013b; McDaniel et al., 2014) and fresh waters (López-Bueno et al., 2009; Roux et al., 2012a), stromatolites (Desnues et al., 2008), confined aquifers (Smith et al., 2013), and pelagic sediments (Yoshida et al., 2013). Based on these observations, degenerate primers designed to amplify fragments of the genes encoding the replication initiator (ORF4 or Rep; Tucker et al., 2011) and MCP (Labonté and Suttle, 2013a; Hopkins et al., 2014) were used to examine the distribution and diversity of gokushoviruses in marine ecosystems. Phylogenetic analyses of the amplified fragments revealed that ssDNA phages have different geographic distributions (Labonté and Suttle, 2013a), and that the genetic distance of gokushovirus sequences increased with geographic distance (Tucker et al., 2011). Most of the hosts of environmental gokushoviruses are unknown, but it is hypothesized that they occupy specialized niches, and that specific gokushovirus genotypes have limited geographic range (Angly et al., 2006; Tucker et al., 2011; Labonté and Suttle, 2013a).