The Arctic virome is dominated by unknown and small ssDNA viruses

On average, blastx searches against known viral genomes could only assign 9.8% of the reads to a previously known taxonomic unit (Fig. 2). The vast majority of these were classified as ssDNA viruses (86%), with an important fraction of Circoviridae (38.1%), followed by dsDNA (double-stranded DNA) viruses (2.8%), mainly composed of Caudovirales (1.8%). Unexpectedly, a fraction of the assigned reads of these DNA viromes were assigned to ssRNA viruses. Comparisons of the taxonomic profiles based on classified reads from the Arctic viromes showed that, overall, they harbored members of the same viral families. However, it is noteworthy that the amount of classified reads was much lower for sample Lv1 than for the other samples (2.5% versus 7.8 to 15.5%) and that the profiles of samples Lv1 and SvL2 separated from the rest mainly because of increased abundances of Circoviridae- and Nanoviridae-classified reads and a lower proportion of Microviridae. The Phi29 polymerase has been shown to preferentially amplify circular ssDNA (17). However, on the basis of the percentage of assigned reads to each of the most abundant viral groups (ssDNA viruses, Siphoviridae, Podoviridae, and Myoviridae) and accounting for the average differences in genome sizes, the relative abundance of ssDNA viruses in the Arctic viromes, assuming the reported 100 × (17) positive amplification bias, was 88.6 ± 8.3%. Therefore, ssDNA viruses were abundant in the Arctic viromes even considering a Phi29 polymerase bias.

Open in a separate window

Fig. 2

Taxonomic distribution of assigned metagenomic reads from Arctic viromes.

Results represent blastx hits against a viral genome database summarized by MEGAN (Meta Genome Analyzer). Only the 10 most abundant categories are shown.

Most polar contigs do not bear high similarity to known viral genomes

Assembly of polar reads, including those with no similarity in databases, into larger sequences produced an average of 15,541 contigs per sample [5175 at least 1000 base pairs (bp) in length], accounting for a total of over 97 Mbp, including contigs as large as 114,603 bp (Table 1). blastx searches revealed that for all Arctic viromes, contigs with the highest abundance were likely ssDNA viruses (fig. S1A). We cross-aligned the contigs obtained to all known viral genomes. No similarities greater than 65% were found for IR1 (Isfjord Radio 1), IR2, and Lv1. Both SvL1 and SvL2 produced contigs (SvL1-4380 and SvL2-1027) with resemblance to Sclerotinia sclerotiorum hypovirulence–associated DNA virus 1 (507-bp overlap, 79% similarity; and 975-bp overlap, 89% similarity, respectively), a fungus virus related to the Geminiviridae. Also, the resequenced Antarctic virome produced a contig (Ant-0) similar to Bathycoccus sp. RCC1105 virus BpV1 (500 bp, 76%), a phycodnavirus, and another contig (Ant-53) with regions similar to Cafeteria roenbergensis virus BV-PW1 (1439 bp, 87%) and Phaeocystis globosa virus (3024 bp, 93%), members of the Mimiviridae and Phycodnaviridae, respectively. Taxonomic affiliation of the contigs using the METAVIR pipeline confirmed the prevalence of ssDNA viruses and the presence of sequences related to ssRNA viruses within the Arctic viromes. A closer inspection of these sequences revealed that they belonged to small circular contigs composed of replication proteins most similar to Circoviridae-Geminiviridae-Nanoviridae genes (and, in a few cases, Phycodnaviridae) and coat proteins most similar to the ssRNA virus Sclerophthora macrospora virus A.

The Arctic freshwater viral communities share genes and genomes

A stringent mapping analysis [minimum 98% similarity, shown to differentiate between closely related phages (18)] of reads to contigs assessing the extent of genetic content overlap between communities (fig. S2) showed very little fine-grain genetic overlap between Arctic and Antarctic viral communities. Lv1 and SvL2 shared the most genetic information. IR1 and IR2 were also genetically very similar to each other, and overall, SvL1 shared the least genetic information with the other Arctic viromes.

Next, we studied to which degree these Arctic freshwater environments might share common species (genomes) and/or highly similar genomic regions, analyzing cross-contig alignments with 98% similarity and 500-bp length thresholds. IR1 and IR2 shared the greatest number of high-similarity genomic regions (table S1), whereas Lv1 and SvL1 shared relatively few of these regions when compared with the rest of the pairwise comparisons. A more detailed analysis detected three groups of highly similar circular contigs (indicative of completion of a genomic element) present in four of the five Arctic communities (fig. S1B). These sets had intra-group similarities >97% and lengths ranging from 1120 to 1710 bp, presenting two to three open reading frames (ORFs). As for the most abundant circular contigs in the Arctic viromes (fig. S1A), blastx searches showed significant similarities for some of these ORFs to known replication proteins of ssDNA viruses, in most cases related to uncultured marine viruses.

The Arctic freshwater virome differs from other environments at fine scale

Analysis of blastx-derived taxonomic profiles showed a strong partitioning of the viromes between polar and nonpolar viromes, mainly driven by the relative abundances of small circular ssDNA viruses and phages (Fig. 3). This clustering includes the Lv1Pond and, significantly, the Antarctic spring virome. On the contrary, Antarctic summer virome segregates from the rest mainly because of increased proportions of Phycodnaviridae-related sequences (table S2). Another feature is that the Arctic Ocean community separates from the Arctic viromes. The Arctic viromes were found to contain mainly genes from the categories of phages, prophages, and transposable elements (not shown), which were also dominant in the temperate Lake Bourget and Antarctic spring viromes. The other prominent feature is the enrichment within all Sahara viromes, Antarctic summer, and Arctic Ocean of reads with similarity to genes coding for functions related to various families of metabolic pathways, including photosynthesis, respiration, and stress response. Nevertheless, we did not proceed with further analysis because the overall percentage of assigned reads was likely too low to provide informative conclusions (average, 1.6 ± 2.9%).

Open in a separate window

Fig. 3

Ordination of individual freshwater viromes according to their taxonomic profiles.

The figure represents the first two axes of a principal components analysis based on the taxonomic affiliation of assigned reads. (A) Projection of the individual viromes on the model. Viromes from the same environment are represented by equal forms. (B) Projection of the viral groups with a higher contribution to the model. The position of each group represents its association with each virome. The analysis indicates a strong partitioning of the viromes between polar and nonpolar viromes, driven mainly by the relative abundance of small circular ssDNA viruses and phages. U., unclassified; d, distance.

Reference-dependent analysis of metagenomic reads returned a small proportion of positive hits. Hence, we also used two reference-independent comparison methods, using all metagenomic reads: crAss, a cross-assembly–based fine-grained analysis, and more coarse-grained cross-tblastx comparisons between viromes. crAss results (Fig. 4) indicated that viromes from the same environment share fine-grained genetic information, with unnoticeable overlap between different environments, including Arctic-Antarctica and Arctic freshwater-seawater. Such clustering by environment is also observable in the cross-tblastx results (fig. S3). However, in this case, some degree of coarse-grain genetic overlap between environments is observable. For instance, Aquaculture and Sahara communities shared some relative overlap (averaging 5.6 ± 1.9%), and the Antarctic showed some overlap with the Arctic (averaging 7.7 ± 4%). Two of the Aquaculture and both temperate lake viromes (especially Bourget) showed some overlap with both the Sahara and Arctic samples, which nevertheless showed little overlap among themselves (0.5 ± 0.64). Finally, the reclaimed water viromes showed little overlap to the other, except for Nursery, which displayed overlap with the Arctic viromes (5.3 ± 3%) except for Lv1 and Lv1Pond. Finally, the Arctic Ocean virome consistently showed very low overlap values (0.3 ± 0.3). We also produced per-environment viromes by pooling all available reads for that particular environment. The results of the cross-tblastx analysis of these pools (Fig. 5) are confirmatory of those obtained with the individual communities with regard to the overall coarse-grain genetic overlap between the different environments. They showed even more clearly how the temperate lake viromes present high relative overlap to all other environments (except Arctic Ocean) and vice versa, likely representing a midpoint between the ssDNA virus–dominated polar communities and the other phage-abundant environments.

Open in a separate window

Fig. 4

Reference-independent cross-assembly (crAss) comparison of individual freshwater viromes.

The dendrogram represents a hierarchical clustering (average-linkage method) of Wootters distance values derived from the cross-assembly of freshwater viromes. Viromes from the same environment are represented by equal forms. The analysis evidences very low fine-grain genetic overlap between viromes from different environments.

Open in a separate window

Fig. 5

Coarse-grain genetic overlap between environments.

The heatplot represents the results of cross-tblastx searches between per-environment pooled viral metagenomes. Values (% of total reads with hits with E value <0.001) are registered within the boxes and color-coded. Horizontal labels represent viromes acting as query, and vertical labels represent viromes acting as reference.

Some ssDNA viral groups exhibit a bipolar distribution

The per-environment pooled reads were also used for contig reconstruction to look for shared genomes and/or genomic regions between the different environments. Overall, the results (table S3) reveal little contig overlap between the different environments, although there is a noticeable trend of increasing overlap with decreasing similarity threshold. Further analysis revealed that two of the overlapping contig pairs between the Arctic and Antarctica represented circular contigs showing 93.8 and 90.8% similarity along their complete sequence (fig. S1C). The first group includes a viral_rep gene most similar to a likely Circoviridae sequence and another ORF bearing low similarity to a phage integrase, whereas the second group presents two ORFs, one of them with similarity to a Circoviridae-related putative viral capsid protein. Also, one of the overlapping contig pairs between the Arctic and the French temperate lakes represented a 1936-bp sequence with 95.8% similarity. This sequence presented two ORFs with no similarities and a third ORF with similarity to a phage structural protein. Strikingly, all overlapping contigs between the Antarctic and Sahara environments at the highest similarity, and the single overlapping contig between the reclaimed water environment and both Antarctic and Sahara environments, produced blastx results pointing to DNA mobilization or replication-related proteins related to Acinetobacter sequences. The longest among the overlapping contigs between Antarctic and Sahara environments represented a circular 5594-bp sequence (Sahara_97 has a 180-bp deletion and then 99.9% similarity with respect to Antarctic_92). Further analysis showed that this contig is formed by a 2878-bp segment with 83% similarity to Acinetobacter baumannii E7 pRAY-v2 plasmid and a second segment with 99% similarity to a hypothetical Lactamase_B superfamily protein from Acinetobacter sp. Thus, this circular contig likely represents an antibiotic resistance–carrying plasmid. The second group (contigs linking the reclaimed water, Sahara, and Antarctic environments) was formed by circular contigs ranging from 1739 to 2088 bp in length and with an average similarity of 74.6%. These sequences bear great similarity to uncultured bacterial and Acinetobacter sp. plasmids and to Sphinx 1.76, a nuclease-resistant circular DNA (containing a replicase ORF related to Microviridae) found to co-purify with infectivity in various transmissible spongiform encephalopathies (19). A recent study (20) has found experimental evidence that Acinetobacter sp. DS002 plasmid, most similar (67%) to the sibling Sphinx 2.36 circular DNA (also transmissible spongiform encephalopathy–related and carrying the putative Microviridae replicase), is a phage. These Acinetobacter-related contigs may represent bacterial DNA contamination that profited from a Phi29 overamplification, but these sequences were unrelated to genomic regions and corresponded to plasmids or known or suspected Acinetobacter phages. These sequences were detected in three different environments and produced in three laboratories using different sequencing platforms. Moreover, they represent distinct parental sequences, likely indicating a natural community origin rather than a unique contaminant. All three protocols included nuclease treatments, indicating that the sequences could correspond to genetic material protected within virions. The sequences were not detected in our recent Arctic viromes, nor did our Phi29 amplification–negative controls produce noticeable mass. Thus, the origin of these DNA sequences is unclear; they likely represent abundant, varied, and resistant environmental contaminants (plasmids), but some could correspond to true viral sequences of global freshwater viromes.

Virome extraction

Viral metagenomic DNA was obtained as described (14) with minor modifications. Ninety liters were filtered through a 30-μm nylon mesh and by 0.45-μm tangential flow filtration (TFF) using a Centramate holder (Pall) to remove bacteria and smaller eukaryotes. Viral fractions were concentrated 100 times by 70-kD TFF. Viral stocks were preserved at −20°C/−80°C before DNA extraction.

DNA extraction and sequencing

Frozen stocks were thawed at 4°C and passed through a 25% sucrose cushion by centrifugation for 16 hours at 60,000g and 4°C. The pellets were resuspended in 10 mM tris (pH 8) and 1 mM EDTA and filtered using a 0.45-μm filter. Viral concentrates were then treated with deoxyribonuclease I (DNase I) (500 U ml⁻¹), nuclease S7 (420 U ml⁻¹), ribonuclease A (RNAse A) (100 μg ml⁻¹), and RNAse H (2 U per reaction) for 30 min at room temperature to remove free nucleic acids. Nuclease reactions were stopped with 12 mM EDTA/2 mM EGTA, and then viral capsids and envelopes were disrupted with SDS (0.5%) and proteinase K (200 μg ml⁻¹) treatment. Viral DNA was extracted with phenol-chloroform, and ethanol was precipitated. The resulting DNA was randomly amplified using Phi29 polymerase and modified random hexamers (GenomiPhi HY, GE Healthcare) for 2.5 hours, according to the manufacturer’s instructions. Finally, the nucleic acids were shotgun-sequenced using an Illumina HiSeq apparatus (Parque Científico de Madrid), resulting in 5.4 million to 7.4 million 2 × 100–bp reads per sample. The extraction protocol used successfully limited possible bacterial contamination, as evidenced by transmission electron microscopy and BLAST searches against an all-inclusive 16S ribosomal RNA gene database. For comparative purposes, metagenomic DNA from the Limnopolar Lake summer sample previously reported (14) was also resequenced with the same device obtaining more than 1.8 million 2 × 75–bp reads. Additionally, DNA derived in the same fashion from Lv1Pond was sequenced with a Roche 454 FLX device obtaining more than 228 thousand reads circa 250 bp in length. Because of its relatively small amount of information, this last data set was only used in the comparisons to other published freshwater viromes.

Metagenomic read analysis

A series of filtering and trimming steps were undertaken to remove low-quality reads and bases using prinseq-lite software (37). Taxonomic assignment of reads was carried out by performing blastx (E score <10⁻³) searches using as reference set the most recent release of viral genomes from the National Center for Biotechnology Information (NCBI) containing 4788 sequences, and then summarizing the results with MEGAN 4 (38). To compare the viromes obtained in this study with other published freshwater viromes, we downloaded their sequences from the relevant sources. These included data sets from an aquaculture facility (Tilapia504, Tilapia608, and Tilapia439; SEED accessions 4440412.3, 4440439.3, and 4440424.3) (32), samples from an Antarctic lake [Antarctic summer and Antarctic spring; SRA (Sequence Read Archive) accession SRA008157] (14), viromes from Sahara desert perennial ponds [ElBerbera, Hamdoum, Ilij, and Molamhar; MG-RAST (Metagenomic Rapid Annotations using Subsystems Technology) ID 4446033.3, 4445715.3, 4445716.3, and 4445718.3] (39), data sets from a reclaimed-waters study (Potal, Effluent, Park, and Nursery; SRA accession SRA008294) (40), samples from French temperate lakes (Pavin and Bourget; SRA accession ERP000339) (41), and one virome from the Arctic Ocean (ArcticOcean; MG-RAST ID 4441621.3) (42). All data sets were subsampled to a common depth (the minimum number of sequences in a virome; 39,351) and then trimmed to 100 bp to homogenize sequence length, thus effectively normalizing sampling effort. Taxonomic and functional assignments were carried out as above but in the latter case using NCBI’s nonredundant database and then using the SEED (43) classification in MEGAN. A cross-tblastx (E score <10⁻³) between these data sets was undertaken to study the degree of putative coarse-grain genetic overlap between these environments. Then, the percentage of reads giving above-threshold results were graphically represented as a heat map using the gplots R package. To further study the relationships between the different freshwater viral communities, all reads per environment were pooled. The genetic overlap was studied as above using a cross-tblastx (using subsamples of 252,867 reads, 70 bp in length). To further study the relationships between the subsampled individual viromes, we followed the recently described reference-independent strategy crAss (44). Briefly, all sequences were combined in a single pool and then assembled using idba_ud (see below) generating cross-contigs. Then, reads from each environment were mapped onto these cross-contigs using Bowtie 2 (--score-min L,0,−0.2) (45). The degree of similarity between these data sets was assessed using the Wootters distance formula on cross-contig mapped-read counts as previously suggested (44) and represented as a cladogram with R. Then, several ecological indices were calculated for the subsampled individual viromes using the PHACCS software (46). To this end, the average genome length was assessed using GAAS (47), and their contig spectra were calculated with Circonspect (48) using the Minimo assembler (49).

Metagenomic assembly and analysis

Arctic and resequenced Antarctic reads were assembled into contigs using idba_ud (50). Assignment of reads to contigs was also done with Bowtie 2. Cross-contig alignments and alignments of contigs to reference genomes were carried out using NUCmer (51) (500-bp minimum overlap, similarity thresholds ranging from 65 to 98%). On the other hand, to assess the viral genomic connectivity between the environments studied, their per-environment pooled reads were assembled into contigs as above and cross-contig comparisons (again using NUCmer) were carried out.

We thank A. Quesada and M. L. Avila for helpful comments on the manuscript. Sample collection of Arctic lakes was done as part of the University Centre in Svalbard Master Course AB-327 on Arctic Microbiology (2010–2012), and we thank M. L. Avila, L. Little, and the master students for their help. Funding: This work was funded by the Spanish Ministry of Economy and Competitiveness grant CTM2011-15091-E/ANT. D.A.d.C. was supported by the Marie Curie International Incoming Fellowship grant PIIF-GA-2012-328287, and A.L.-B. was supported by a Ramón y Cajal Fellowship (RYC-2010-06300). Author contributions: D.A.P. and A.A. conceived the project and obtained and concentrated all virome samples. D.A.d.C. and A.L.-B. extracted the metagenomic DNA and prepared the libraries. D.A.d.C. analyzed the sequence data. D.A.d.C., A.L.-B., and A.A. wrote the initial manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: Metagenomic data have been deposited at SRA under accession PRJEB5265, and polar contigs have been deposited at METAVIR under IDs 3295 to 3300.