Dynamics of Viral Abundance and Diversity in a Sphagnum-Dominated Peatland: Temporal Fluctuations Prevail Over Habitat

Viruses impact microbial activity and carbon cycling in various environments, but their diversity and ecological importance in Sphagnum-peatlands are unknown. Abundances of viral particles and prokaryotes were monitored bi-monthly at a fen and a bog at two different layers of the peat surface. Viral particle abundance ranged from 1.7 x 106 to 5.6 x 108 particles mL-1, and did not differ between fen and bog but showed seasonal fluctuations. These fluctuations were positively correlated with prokaryote abundance and dissolved organic carbon, and negatively correlated with water-table height and dissolved oxygen. Using shotgun metagenomics we observed a shift in viral diversity between winter/spring and summer/autumn, indicating a seasonal succession of viral communities, mainly driven by weather-related environmental changes. Based on the seasonal asynchrony between viral and microbial diversity, we hypothesize a seasonal shift in the active microbial communities associated with a shift from lysogenic to lytic lifestyles. Our results suggest that temporal variations of environmental conditions rather than current habitat differences control the dynamics of virus-host interactions in Sphagnum-dominated peatlands.


INTRODUCTION
Viruses are present in virtually all ecosystems, and are considered to be the most abundant biological entities in the biosphere (Suttle, 2005). Through lytic and lysogenic life cycles, viruses affect the metabolism and abundance of their cellular hosts from all three domains of life (synthesis in) (Kirchman, 2012), impacting the diversity and the structure of microbial communities (Middelboe, 2000;Suttle, 2007;Sime-Ngando, 2014). Because viruses function as a top-down control on microbial production (Pradeep Ram et al., 2011;Chow et al., 2014), they affect biogeochemical cycles through the release of substantial amounts of organic carbon and nutrients in the environment (Middelboe and Lyck, 2002). While resources released by this viral shunt can be reused to sustain microbial biomass, it typically results in a reduction of the overall flow of organic matter and energy towards higher trophic levels (Fuhrman, 1999;Middelboe and Lyck, 2002;Ankrah et al., 2014). Besides the viral-shunt, a significant proportion of the released viruses represent a consequential amount of labile organic matter that can be readily decomposed (Dell'Anno et al., 2015). Viral lysis, along with bacterial grazing are important processes in microbial succession (Fierer et al., 2010), suggesting that coupled viral-host interactions can influence ecosystem-level carbon cycling, depending on the activity of cells and on the balance between lytic and lysogenic strategies. Consequently, considering viral ecology is critical to understanding ecosystem functioning.
Peatlands are a globally relevant component of the carbon cycle, storing a quarter of global soil carbon and more carbon than all vegetation (Yu, 2012;Turetsky et al., 2015). Boreal and temperate peatlands are often dominated by peat forming mosses, belonging to the Sphagnum genus. Many Sphagnumdominated peatlands develop from a nutrient-rich, non-acidic fen (minerotrophic fen) during early stages, to a nutrient-poor acidic bog (rain-fed ombrotrophic bog) (Rydin and Jeglum, 2006). During this development process, peatlands accumulate large stocks of partly decomposed organic matter as peat. This accumulation is a consequence of the long-term imbalance between carbon uptake from photosynthesis and carbon losses via respiration, methanogenesis and DOC leaching, which results in the preservation of up to 15 % of the net primary production as peat (Clymo, 1984;Moore, 1987;Francez and Vasander, 1995;Roulet et al., 2007). The combination of abiotic conditions, such as low temperature, low pH, high soil water content and poor nutrient availability, as well as biotic factors (Sphagnum litter quality, allelopathy) constrains microbial activities, limiting decomposition and mineralization (Clymo, 1984;Rydin and Jeglum, 2006). Accumulation of organic matter over long time periods leads to stratification with a deep permanently anoxic catotelm in contact with underlying bedrock, an intermediate mesotelm layer that is periodically anoxic and oxic and the predominantly oxic acrotelm at the surface (Clymo and Bryant, 2008). The heterogeneity of peat, the variety of dynamic stages and sources of water and nutrients provide a large panel of habitats and niches for a broad diversity of microorganisms transforming different carbon substrates through aerobic and anaerobic pathways of decomposition (Juottonen et al., 2005;Artz, 2009;Andersen et al., 2010;Tveit et al., 2012).
Understanding of microbial diversity in peatlands, especially in regards to spatio-temporal patterns, is currently limited. The development of cultivation-independent approaches such as metagenomics provide a powerful tool to investigate microbial taxa and associated protein-coding genes across ecosystems (Lynch et al., 2004;Vandenkoornhuyse et al., 2010). More recently, sequencing of viral metagenomes (hereafter metaviromes) has dramatically expanded understanding of viral diversity, particularly in marine and freshwater environments (Comeau et al., 2008;Lopez-Bueno et al., 2009;Rohwer and Thurber, 2009;Roux et al., 2012;Smith et al., 2013;Labonté and Suttle, 2013;Wommack et al., 2015) compared to soil ecosystems (Kimura et al., 2008;Zablocki et al., 2014;Adriaenssens et al., 2015). Most analyses have focused on cross-biomes comparisons, while spatio-temporal dynamics remain poorly documented (Pagarete et al., 2013;Chow et al., 2014). Several viral marker genes, conserved within particular viral families, are considered to be reliable for taxonomic affiliation (Roux et al., 2011;Sakowski et al., 2014). Nevertheless, none of these marker genes are ubiquitous among viruses and they are far from the high level of conservation found in rRNA genes that are used for the classification of organisms from the three kingdoms.
While several studies have investigated the diversity of microbial communities in peatlands (Dedysh et al., 2006;Artz et al., 2007;Peltoniemi et al., 2009;Dedysh, 2011;Mackelprang et al., 2011;Bragina et al., 2012;Serkebaeva et al., 2013;Mondav et al., 2014;Tveit et al., 2014;Hultman et al., 2015;Nunes et al., 2015) viral diversity in Sphagnum-dominated peatlands remains largely unknown and basic knowledge of virus ecology in these ecosystems is still lacking . In view of the importance of viruses in structuring and regulating prokaryotic communities and the implication of the latter in the carbon sink function of peatland; it is essential to understand the role of viruses in the dynamics of the microbial communities in this ecosystem.
In order to characterize virus ecology in a temperate Sphagnum-dominated peatland, we combined and integrated different approaches to study the spatio-temporal patterns of viral abundance and diversity. The goals of this work were: (1) to compare the seasonal abundance of viruses and prokaryotes at two different layers of the peat surface corresponding to a stratified analysis of acrotelm, the most active layer, of fen and bog; (2) to identify abiotic controls on abundance and distribution of these viruses and prokaryotes; and (3) to characterize the viral diversity over an annual cycle.

Sampling Strategy and Experimental Design
Twelve field sessions were organized between May 2010 and November 2012 (Supplementary Table S1): One in May 2010, four in 2011 (June, August 12th, August 24th and October), and seven in 2012 (March, May, June, July, August, September, and November).
For each field session, three peat profiles (replicates) were analyzed in the fen and the bog. In order to avoid disturbance during the extraction of peat and water, sampling for viral counts and metagenomes followed a 3-step progressive cutting from surface to depth. The 0−5 cm layer corresponding to the living Sphagnum carpet (capitulum and "green" stems and leaves) was cut off and carefully removed. The 5−10 cm layer (called "upper surface layer" thereafter) was then cut off and a fraction of peat matrix was harvested and stored in a 50 mL tube for metagenome production. From the remaining undisturbed peat, pore-water was expressed for viral and prokaryote counts. The same method was applied to the 10−15 cm layer (called "lower surface layer" thereafter). For metaviromes production as well as for physicochemical analysis pore-water was expressed from the remaining peat surface layers combined. At least 200 mL of pore-water were necessary for metavirome production.

Physico-Chemical Parameters
Peat temperature was measured directly in the peat profiles at 5, 10, and 15 cm under the Sphagnum capitula layer. Conductivity, pH and oxygen were measured in the field with filtered (125 µm) peat water. Dissolved organic carbon (DOC) and anion concentrations (nitrate, sulfate, and acetate) were measured at the "Biogeochemical Analysis" platform (ECOBIO and GEOSCIENCES -OSU Rennes), following water filtration at 0.2 µm (Whatman), using a Bioritech DOC Analyser and a Dionex Analyzer (Table 1).

Metavirome and Metagenome Production
To get enough material fen and bog viral communities were sampled from pore water combining the upper and lower surface layers. Samples were recovered at four dates (Supplementary Table S1, Table 2): 07 June 2011 (vFen_June11, vBog_June11; along with peat for a study of the microbial communities), 12 August 2011 (vFen_Aug11, vBog_Aug11), 12 October 2011 (vFen_Oct11, vBog_Oct11), and 23 March 2012 (vFen_Mar12, vBog_Mar12). Due to low water content of the peat, only one sample was collected for each dynamic stage for the three first dates. In March 2012, the fen and bog were sampled in biological triplicates (vFen_Mar12_A, B, C; vBog_Mar12_A, B, C). In total we prepared six fen and six bog metaviromes. Sphagnum water was prefiltered at 125 µm. Viruses were concentrated using PEGylation (Colombet et al., 2007). Viral concentrates were filtered through a 0.20 µm filter (Minisart, Sartorius) and diluted 5x in H 2 O (Sigma) to a volume of 5 mL. Remaining nonviral DNA was digested with 10 U DNAse RQ1 RNAse free (Promega) at 37 • C for 1 h. Viral DNA was extracted as described previously . DNA quality was checked with the High Sensitivity DNA kit on a 2100 Bioanalyzer (Agilent). Whole genome amplification (WGA) was applied in triplicate for each sample using GenomiPhi V2 (GE Healthcare) following manufacturers' instructions and the triplicates were pooled.
Subsequent pyrosequencing was performed on a GS FLX system (454 Life Sciences, Roche, Branford, CT, USA) at the "Functional and Environmental Genomics" platform (OSU, Rennes, France). Roche/454 filtering tools and stringent filters were used to ensure the highest sequence quality and to remove artificial replicates of sequences and sequences smaller than 250 bp as shown previously (Quaiser et al., 2014).
For the metagenomes, 12 un-pressed peat matrix samples, distinguishing bog/fen, upper and lower surface layers, each in triplicates, collected on 07 June 2011, underwent DNA extraction, pyrosequencing as well as size and quality trimming (Supplementary Table S2) as described previously (Quaiser et al., 2014). Briefly, 2 g of peat matrix were lysed in 15ml of extraction buffer containing 4% cetyltrimethylammonium bromide (CTAB), 0.5% polyvinylpyrrolidone (PVP, Sigma−Aldrich), 0.7 M NaCl, 100 mM potassium phosphate (pH 6.8), 20 mM EDTA (pH 8.0), 1% beta-mercaptoethanol, 1 M guanidin thiocyanate and incubated at 65 • C for 30 min. Homogenization was obtained by vigorous vortexing every 5 min during 1 min in the presence of glass beads. One volume of chloroform-isoamylalcohol (24:1) was added, vortexed for   1 min and incubated at room temperature for 5 min. The samples were centrifuged at 4000 g for 15 min at 4 • C and the aqueous phases were transferred to new tubes. Binding conditions for silica-based RNA extraction were adjusted, applied on Nucleo Spin RNA II kit columns and subsequent purification was performed following the instructions of the manufacturer (Macherey−Nagel). DNA was nebulized to fragments of about 700 bp. The DNA was purified with Agencourt AMPur XP magnetic beads (Beckman−Coulter). DNA fragmentation quality was checked with the High sensitivity DNA kit on a 2100 Bioanalyzer (Agilent). Subsequent library construction and pyrosequencing was performed in technical duplicates on a GS FLX system (454/Roche) at the "Functional and Environmental Genomics" platform (OSU, Rennes, France). Roche/454 filtering tools and stringent filters developed locally were used to ensure the highest possible sequence quality and to suppress artificial replicates of sequences as well as sequences smaller than 250 bp. The metaviromes are available under the Metavir IDs (http://metavir-meb.univ-bpclermont.fr/): 1368, 1369, 1370, 1373, 1374, 1375, 1376, 1377, 1378, 1379, 1380, and 1382. In addition, pyrosequencing reads reported in this publication have been deposited in the ENA Sequence Read Archive under the study accession number PRJEB11420 (metaviromes) and PRJEB11421 (metagenomes).

Metavirome Analysis
After pyrosequencing, sequence quality and size trimming, we obtained 481 402 and 615 487 sequences with an average length of 415 bp from fen and bog, respectively ( Table 2). The quality of the virome extraction process was assessed by determining the amount of rRNA and tRNA sequences using Meta_RNA that identifies SSU and LSU rRNAs from the three kingdoms (Huang et al., 2009) and tRNAscan-SE (Lowe and Eddy, 1997). In total 69 rRNA (0.0063%) and 1681 tRNA (0.15%) sequences were identified indicating an insignificant level of potential contamination of microbial DNA. To avoid misinterpretation of the results, these sequences were excluded in subsequent analysis.
The viral diversity was analyzed using Metavir (Roux et al., 2014) and the sequences were subjected to tBLASTx (Altschul et al., 1997) against the NCBI RefSeqVirus database (e-value 10 −7 ). Taxonomic assignment of the sequences was determined with MEGAN (Huson et al., 2007). Several accompanying tools were used on the Galaxy/Genouest bioinformatics platform (Le Bras et al., 2013).
To estimate the level of similarity between the viral communities, the proportion of similar sequences of each pair of metaviromes was computed with Compareads 1.2.2 (Maillet et al., 2012). Using this software, two sequences are considered to be similar if they present a defined number (t) of identical k-mers (k). To calibrate this analysis we tested 3 different numbers (t = 2; 4; 10) of identical 33mers (k = 33 nt). The most reliable results were obtained using four identical 33mers, parameters that were used for further analyses computed with Compareads. Compareads output is a percentage of similarity between a pair of metaviromes. These percentages were used to build a distance matrix, on which hierarchical clustering was performed using the R package pvclust package (Suzuki and Shimodaira, 2006) (distance = "correlation", method = "average").
Comparisons made with Compareads give a global estimate of the similarity between the metaviromes. In order to take into account the diversity of sequences within each metavirome, we analyzed the qualitative distribution (presence−absence) of clusters of highly similar sequences in the 12 metaviromes. This second analysis also allowed removal of potential bias due to the variation of the number of sequences obtained for each metavirome. Sequences from the 12 metaviromes were clustered using CD-HIT-EST (Huang et al., 2010) (c = 0.95; n = 8). Clustering results were used to compute Sørensen dissimilarity between pairs of metaviromes using MOTHUR (Schloss et al., 2009). Hierarchical clustering (pv-clust package) was used to represent compositional relatedness between metaviromes from the matrix of Sørensen dissimilarities. Clusters of sequences were split into different categories according to the amount of sequences they contained and the same analysis was performed for each size category.

Analysis of the Sequences Shared by a Metagenome and a Metavirome
In order to find viral sequences in the metagenomes, and to analyze the link between the viral and microbial communities, sequences shared by at least one metagenome and one metavirome were retrieved using Compareads (k = 33, t = 4), and clustered using CD-HIT-EST (see "Metaviromes Analysis"). To ensure that this selection of sequences did not alter the compositional patterns observed for the total metaviromes, Sørensen dissimilarity was calculated with these clusters for the "shared" sequences originating from the metaviromes. Correspondence Analysis (CA) was performed on the whole dataset of shared sequences and the sample dissimilarities carried on the two first axes was represented with a hierarchical clustering (ade4 package) (Dray and Dufour, 2007). Taxonomic assignment of the sequences was obtained using Metavir tBLASTx output (e-value 10 −7 ) (Roux et al., 2014), and analyzed with MEGAN (Huson et al., 2007).

Statistical Analysis
Due to non-homogeneity of variance, one-factor Kruskal-Wallis tests were used on PA, VPA, and VPR in order to detect differences between sites (fen vs. bog), layers (upper vs. lower) and sampling dates. A principal component analysis (PCA) was performed on the physico-chemical dataset taking into account samples with available DOC (ade4 package) (Dray and Dufour, 2007). The first component was associated with the fluctuation of the physico-chemical variables through the habitats and seasons. We used the sample coordinates on the first component as a variable representing the spatio-temporal gradient (Legendre and Legendre, 1998;Ramette, 2007). Then potential relations between the gradient and log transformed PA, VPA, and VPR were tested using linear regression. All statistical analyses were performed using the open-source statistical software R (version 2.14) (R Development Core Team, 2013).

Link Between Abiotic Variables, VPA and PA
The potential relationships between abiotic and biotic environment and virus communities were analyzed combining fluctuation of physico-chemical variables with viral and prokaryote abundance. The annual mean water-table level was −2.9 ± 4.5 cm in the fen and -18.0 ± 5.0 cm in the bog. Temperature varied more at the upper layer of the bog (annual mean: 11.6 ± 4.4 • C) than at the upper layer of the fen (12.4 ± 3.8 • C), but was lower at the lower layer of both stages (fen: 9.3 ± 3.5 • C; bog: 9.2 ± 3.5 • C) ( Table 1). We characterized the variation of abiotic physico-chemical parameters with a (PCA) (Figure 3). The first component of the PCA accounted for 49.8% of the variance, and was positively correlated with conductivity, sulfate (SO 4 2− ) and DOC and negatively correlated with water-table level and oxygen (O 2 ). The second component accounted for 16.4% of the variance and was mainly correlated with nitrate concentration, which was higher in the fen in June, July, and September. The distribution of data points across the first two components emphasized the differences between fen and bog but also highlighted similar temporal trends within the two development stages, distinguishing the March and May samples from September. The seasonal fluctuations of the water-table were closely linked with variation in water-chemistry, suggesting that water sources and flowpaths affect nutrient concentrations, potentially due to dilution of a limited solute stock (for example, Spearman's rank correlations; water-table to conductivity: R = −0.67, P = 6 × 10 −5 , N = 29; water-table to sulfate: R = −0.66, P = 9 × 10 −5 , N = 29) (Table 1; Figure 3). Because the first component represents an integrated variable of spatiotemporal variations of the physico-chemical and hydrological parameters, we used sample coordinates on this axis to model the seasonal abiotic fluctuations for both fen and bog. Log-transformed PA and VPA were both positively, linearly correlated with this abiotic gradient (N = 29, PA: r 2 = 0.53, P < 10 −5 ; VPA: r 2 = 0.41, P < 10 −3 ; Supplementary Figure  S1), whereas VPR was unrelated to this gradient (r 2 = −0.03, P = 0.68).

Viral Community Composition and Diversity
Pyrosequencing of the 12 metaviromes yielded 1 096 889 sequences with an average length of 415 bp ( Table 2). Analysis revealed that sequences of ribosomal RNA genes accounted for 0.0063% (69 rRNA sequences), which were excluded from subsequent analysis. In addition, the predicted protein coding sequences that matched the functional category databases (KEGG) and the cluster of orthologous genes database (COG) accounted for only 0.38 and 0.45%, respectively. This corresponds to 10−50 times fewer matches than are typically found in conventional short read metagenomes ( Table 2; Quaiser et al., 2011), indicating a very low level of contamination by genomic DNA from microorganisms (Roux et al., 2013b). This allows the precise characterization of viral diversity and variation in fen and bog through the year.

Taxonomic Composition
Sequences were compared against viral genomes from the NCBI RefSeqVirus database. Only a small proportion of sequences, ranging from 4.2% (vBog_Oct11, v = virus/Fen or Bog/sampling date) to 10.9% (vBog_June11) matched the available viral genomes indicating the presence of currently undetected viruses (Table 2; Figure 4). Matches associated with ssDNA viruses were most common, accounting for a mean of 4.5% of the total number of sequences, with primary assignment to the bacteriophage family Microviridae (1.7 ± 1.3%) and to the eukaryal ssDNA family Circoviridae (0.9 ± 0.7%). Matches with dsDNA viruses appeared mostly affiliated with the order of Caudovirales (1.3 ± 0.8%), which can be hosted by both Bacteria and Archaea. While the protocol was not designed to preserve RNA viruses, we detected sequences matching to ssRNA viruses affiliated with Tombusviridae and Sclerophtora macrospora virus A representing likely the recently identified so called "chimeric viruses" (Diemer and Stedman, 2012;Roux et al., 2013a). They were present in all samples and accounted for 0.1% (vFen_Aug11) to up to 2.8% (vBog_Aug11) of the total metavirome sequences (Figure 4). Nevertheless, the interpretation of these results must be considered with respect to the applied multiple displacement amplification, that was shown to be bias prone towards ssDNA viruses (Kim and Bae, 2011). Due to potential biases no statistical analyses were performed on the proportions of the viral types.

Genomic Diversity Based Comparisons of the Metaviromes
Due to the lack of virus reference sequences in the databases the majority of sequences (88.3 to 95.4%) remained unaffiliated to known viral taxa. To characterize the remaining, unidentified metaviromes sequences, we analyzed the proportion of similar sequences (four identical kmers of 33 bp) between each pair of metaviromes with Compareads (Maillet et al., 2012). The Sequence comparisons with Compareads provide a global estimate of the proportion of similar sequences without taking into account the internal structure of the sequence sets. Therefore we clustered sequences with a 95% identity threshold to assess the diversity of protein-coding gene sequences and to determine which groups of sequences drive the similarity between metaviromes in the Compareads analysis. Sørensen dissimilarity was calculated for every pair of metaviromes. This distance is based on the distribution (presence-absence) of sequences Frontiers in Microbiology | www.frontiersin.org from the different clusters in the 12 metaviromes. Sørensen dissimilarities between pairwise metaviromes were uncorrelated with the amount of sequences in pairwise comparisons (Spearman's rank correlation, R = 0.08, P = 0.495, N = 66). The pattern of similarity between the metaviromes supported the first analysis done with Compareads (Supplementary Figure  S2B), with nearly identical grouping of summer and autumn communities and winter and spring communities.
Unique sequences represented 28% of the total metaviromes and 49% of the sequences belonged to clusters consisting of 10 to more than 1 000 sequences. We split clusters of highly similar sequences into different size categories, depending on the number of sequences they included ( Figure 5) and computed hierarchical clustering based on Sørensen dissimilarities for each size category. Dendrograms indicated the same contrasted pattern between March viral communities and the summer and autumn group of metaviromes to the exception of vFen_Aug11. This distinction was not significant for clusters smaller than five sequences. This suggests a fundamental change in viral community between the two main groups of metaviromes with low intergroup Sørensen similarity (Supplementary Figure S3). Resemblance between winter and spring metaviromes (Fen and Bog from June 2011 and March 2012) was only significant for clusters larger than 250 sequences. Thus, the resemblance between June and March metaviromes appears to be due to a small number of large clusters.
To analyze the genetic similarities with other metaviromes, we compared the 12 peatland metaviromes with 49 available metaviromes from eight different ecosystem types by hierarchical clustering and tBLASTx (Supplementary Figure S4, see Materials and Methods). Peatland viruses formed a distinct group, clearly separated even from geographically close viral communities originating from freshwater lakes, indicating that these metaviromes represent a unique community characteristic of and structured by its ecosystem (Roux et al., 2012).

Link Between Viral and Microbial Communities
To investigate the interactions between viruses and the microbial communities in Sphagnum-dominated peatlands, we sequenced 12 metagenomes from the fen and bog prokaryotic communities (Supplementary Table S2) from the same day and site as the metaviromes vFen_June11 and vBog_June11 (Supplementary Table S1). Metagenome DNA was extracted from the peat matrix allowing finer spatial sampling. In addition, the peat matrix contained the peat pore-water from which the viral particles were sampled. Based on taxonomic affiliations fen and bog prokaryotic communities (hereafter called pFen and pBog) appeared to be predominantly composed of the same main phyla (Supplementary Figure S5). However, regardless of depth, fen and bog appeared to harbor prokaryotic communities with distinct structures as shown by non-metric multidimensional scaling ordination of euclidean distances between metagenomes (Supplementary Figure S6) and with the analysis of similarity (ANOSIM, Euclidean distance, R = 0.68, P < 0.01). In order to identify viral signatures in these metagenomes, we identified sequences shared by metagenomes and metaviromes using Compareads (four identical kmers of 33 bp). We obtained 18,676 "shared" sequences, of which 30 were from bog metagenomes (pBog), 320 were from fen metagenomes (pFen) ( Supplementary  Table S2), and 18,326 from the 12 metaviromes (Table 2). In most metaviromes, the number of "shared" sequences represented less than 1% of the total sequences but reached up to 12% in vFen_Oct11 ( Table 2). We clustered sequences with a 95% identity threshold, built a contingency matrix based on the number of sequences from each sample in the different clusters and performed a CA on this contingency matrix. The first two axes of the CA contained 21% of the total information. Hierarchical clustering based on the first two components revealed three groups of samples (Figure 6). Two groups included metagenomes sampled in June 2011 from fen and bog, respectively (pFen and pBog_ samples) as well as metaviromes from August and October 2011. The third group contained the metaviromes from June 2011 and March 2012. Thus metaviromes from June 2011 did not cluster with the metagenomes sampled on the same day. Clustering and Sørensen dissimilarity analysis based solely on the subset of metavirome sequences shared with metagenomes revealed the same summer and autumn and winter and spring groupings as obtained for complete metaviromes (Supplementary Figure S7). Among these 18 676 "shared" sequences, a total of 774 sequences (4%) were assigned to references in RefSeqVirus (tBLASTx, e-value 10 −7 ) indicating that the vast majority (96%) originate from currently unidentified viruses. Most hits were associated with the Microviridae subfamily Gokushovirinae (ssDNA viruses, 409 hits) and to a lesser extend with the ssDNA viruses Circoviridae (173 hits), Caudovirales (dsDNA viruses, 134 hits), and Sclerophtora macrospora virus A (ssRNA viruses, 19 hits).

Quality of the Metaviromes
To explore the diversity and potential ecological role of viruses in Sphagnum-dominated peatland, we analyzed and compared six fen and six bog metaviromes covering the seasonal periods. The quality of the viromes is essential for comprehensive analysis, since contaminations with microbial genomic DNA would falsify the results. To assure that microbial DNA contamination was satisfactorily low, we applied the pegylation procedure to enrich viral particles (Colombet et al., 2007), DNAse treatment to degrade "free" DNA not protected by capsids, and triplicate whole genome amplification to balance potential amplification bias. The high quality of the virome sequences was shown by the very low abundance of rRNA sequences as well as the low number of matches to functional databases (i.e., KEGG_KO and COG). This is in accordance with the high diversity present in the viral genomic pool and with the high rates of evolutionary changes in viral genomes that are much less conserved than microbial genes (Duffy et al., 2008). Whole genome amplification, as applied here, is known to amplify preferentially circular ssDNA viruses (Kim and Bae, 2011), a bias that cannot be prevented when pooling separate triplicate amplification (Marine et al., 2014). Therefore, the interpretation of the results must be considered with caution. Nevertheless, since all viromes were generated the same way, the inevitably introduced biases should be the same for all allowing reliable comparative analysis.
It has been hypothesized that viruses infecting eukaryotes might be more important in terrestrial ecosystems and wetlands, where protozoan and fungal biomass is higher, while bacteriophages dominate viral consortia in marine and freshwater ecosystems (Farnell-Jackson and Ward, 2003;Jackson and Jackson, 2008;Kimura et al., 2008). Concerning the peatland metaviromes, among the sequences that matched viral genome databases, similar proportions of sequences were assigned to eukaryotic viruses, such as Circoviridae and Sclerophtora macrospora virus A-like viruses, and to prokaryotic viruses, such as Caudovirales or Microviridae. However, due to the vast majority of sequences being unassigned and to the high variability of viral genes, it remains impossible to determine whether viral communities in peatlands are dominated by prokaryote or eukaryote infecting viruses.

Successional Patterns of Viral and Microbial Communities
Ecological integration of the viral compartment into ecosystem functioning is mainly obtained through approaches combining virome sequences and viral abundance analysis (Wommack et al., 2015). In order to characterize viral ecology of Sphagnumdominated peatlands, we monitored seasonal abundance and diversity of viruses and prokaryotes at two different depths of fen and bog and attempted to identify whether these were correlated to abiotic factors. Fens and bogs are development stages of peatlands that differ fundamentally in vegetation (Sphagnum and vascular plants) and associated physico-chemistry (Rydin and Jeglum, 2006). Thus, as already observed for microbial diversity (Opelt et al., 2007;Bragina et al., 2012), we hypothesized that peatland development stage and associated Sphagnum habitat would be the major driver in the distribution of virus and prokaryote abundance and diversity. Our results confirm that the physico-chemical conditions and the structure of the prokaryotic communities differ between the 2 dynamic stages, but, surprisingly, we did not observe any significant difference in VPA, PA, and viral diversity between fen and bog. The viral communities showed no systematic spatial trend and high variability even within replicates.
While we did not detect any significant spatial differences, we observed a significant seasonal fluctuation of virus diversity and abundance. For both fen and bog, VPA and PA (log-transformed data) were strongly correlated with the seasonal fluctuations of water-table, DOC, conductivity and sulfate: i.e., VPA and PA were higher when water-table was low and DOC and sulfate were high (Supplementary Figure S1). DOC has been recognized as a key factor in the C-balance of Sphagnum-dominated peatlands (Billett et al., 2004) and its patterns are driven by both biological activity (microbial production and consumption, plant exudation) and abiotic variables such temperature, water-table level or acidity (Clark et al., 2009) with seasonal fluctuations as a consequence (Moore, 1987). While the pH is recognized as an integrated physico-chemical variable, we did not register strong influence on PA and VPA. The increased pH observed during June and July in the fen is potentially due to photosynthetic activity with strong assimilation of dissolved inorganic carbon by microalgae, which significantly develop in submerged Sphagnum-fen at the beginning of summer (Gilbert et al., 1998). As viral activity is dependent on bacterial production (Middelboe, 2000), the correlations between DOC and VPA and PA suggest a net production of DOC with increasing microbial activity including viral lysis. DOC concentrations also depend on temperature and water flows, which play a key role in the production and redistribution of carbon in the peat (Waddington and Roulet, 1997;Clark et al., 2009), the two factors interacting during drawdown and flooding periods. However, we did not evidence a clear combination of DOC and temperature. The fluctuations of the water-table depend largely on hydrologic inputs that are rainfall in the bog and mainly runoff and groundwater inflow in the fen, which occur at the yearly scale (seasons) and at a short-time scale (episodic events). Sulfate also interacts with DOC concentrations, especially during drought periods when significant production of sulfate may increase peat acidification and ionic strength (Clark et al., 2009). In our study, sulfate concentrations were significantly higher during watertable drawdowns and were positively linked to DOC. In addition, we could not show whether the production of sulfate was linked to nitrate through microbial sulfur oxidation as already demonstrated in peatlands (Burgin and Hamilton, 2008). Nitrate seems to play a key role only during summer in the fen (see component 2 of the ACP), a dynamic stage in which nitrogen mineralization significantly occurs at this season (Francez and Loiseau, 1999). Our results suggest that viruses and prokaryotes are more abundant at the lower surface layer regardless of peatland development stage. This is in accordance with a previous study concerning prokaryote abundance (Dedysh et al., 2006) and potentially due to more buffered temperature and water-table fluctuations at the lower layer providing more stable conditions. Altogether, these results indicate that seasonal changes in temperature and precipitations (allogenic variables) influence PA and VPA via water-table fluctuations and consequently nutrient concentrations with larger effects at the upper than at the lower layer.
Comparisons of metaviromes with metagenomes showed that temporal variations are more influential than differences in peatland habitat in structuring viral communities. There was a particularly substantial shift in sequence composition from spring to autumn with distinct patterns of composition and abundance, suggesting the existence of ecological succession of viral communities at the seasonal scale. This pattern appears to be consistent inter-annually, with metaviromes sampled about a year apart (June 2011 and March 2012) clustering together. These findings suggest, for the first time, that a cyclic succession in peatlands affects free-occurring viruses at the community level. Recent studies on marine ecosystems also described seasonal fluctuations of viral communities at the ocean surface (Chow and Fuhrman, 2012;Pagarete et al., 2013). In these studies, which focused on the diversity of a viral gene marker, seasonality was mainly characterized by fluctuations of dominant viral types while in our study, seasonality was associated with a general change in the composition of viral communities.

Viruses and Carbon Cycling in Sphagnum-Dominated Peatlands
Viruses are believed to be key components of the carbon cycle in many ecosystems, both altering carbon fluxes and contributing to C-redistribution through bacterial lysis (Fuhrman, 1999;Middelboe and Lyck, 2002;Ankrah et al., 2014). Despite recent analyses of peatland microbial food-webs (Lamentowicz et al., 2013), the significance of viruses in the functioning of Sphagnumdominated peatlands remains unknown.
In viral ecology, VPR is generally considered as an indicator of the bacterial hosts metabolic state (Williamson et al., 2005;Kimura et al., 2008) because viral burst size, and thus viral abundance is positively correlated with microbial growth rate (Middelboe, 2000). In the studied peatland, the VPR was low and did not differ between fen and bog, despite the differences between dead organic matter produced in the two dynamic stages, in relation to the dominant Sphagnum species (Francez, 1995;Thormann et al., 2003). The low VPR compared with other ecosystems is likely due to lower metabolic activity of microorganisms, that is in accordance with the functioning of Sphagnum-dominated peatlands where decomposition is slowed down due to constraining conditions (Rydin and Jeglum, 2006;Artz, 2009) and the presence of a significant proportion of dormant cells in the community (Dedysh et al., 2006;Pankratov et al., 2011).
We detected a VPR peak in June 2011, just before a broad modification of viral community composition in fen and bog. Viruses interact with their hosts through at least two main strategies: the lytic and the lysogenic life cycles, the latter is believed to be favored when microbial activity is low (Danovaro et al., 2002;Payet and Suttle, 2013;Sime-Ngando, 2014). The change in the viral community composition in summer could result from a seasonal shift in the active part of the microbial community and related C-cycling processes via decomposition that show seasonal patterns (Basiliko et al., 2005;Sun et al., 2012). This illustrates a transition from lysogenic to lytic strategies of the viruses infecting the newly active prokaryotes. This hypothesis is supported by the low PA associated with the VPR peak in June 2011, which could result from virus-mediated bacterial lysis, and by the similarities between spring metagenomes (June) and summer and autumn viromes, suggesting the presence of prophages in the microbial genomes in June, that were later released and detected in the metaviromes in August and October.

CONCLUSION
We applied an integrated approach linking virome sequence analysis, viral particle and prokaryote abundance, physicochemical parameters and metagenome-virome comparison to get insights into the ecological functioning of the viral community in peatlands. We found that viral community abundance and diversity in Sphagnum-dominated peatlands express an ecological succession, and that viruses, as well as their hosts, are strongly influenced by the temporal fluctuations at the peatland surface. The observed low VPR, compared with other ecosystems, is in accordance with slowed down decomposition processes in Sphagnum-dominated peatlands. The observed shift in viral diversity suggests a seasonal change of the microbial community and the associated switch of viral life-cycle strategies during summer and autumn highlighting the importance of virushost interactions as they control the dynamics of microbial communities. These patterns may be related to changes in C-cycling processes but further studies are needed to strictly link microbial and viral diversity with C-transformations in the peat. These should focus on the seasonality of viruses that infect different kind of hosts (prokaryotes/eukaryotes) in order to identify the main factors driving this succession and changes in functioning processes, as already suggested for seasonal fluctuations of plankton (Sommer et al., 2012).
To assess the ecological relevance of the diversity of viral communities, analysis of the spatio-temporal dynamics of ecosystem-specific metaviromes as applied here, rather than cross-biomes comparisons, represent a powerful approach to overcome the lack of viral genomes in the databases, and to take advantage of the whole diversity carried in the sequenced viral communities.

FUNDING
This work was supported by funding from the INSU-CNRS EC2CO program to A-JF ("Virtou") (Program MICROBIEN, 2011 and grants from the University of Rennes/CNRS and from the "Observatoire des Sciences de l'Univers de Rennes"(OSUR) to AQ.