Viromes outperform total metagenomes in revealing the spatiotemporal patterns of agricultural soil viral communities

Viromes outperform total metagenomes in the recovery of viral sequences from complex soil communities

To determine the extent to which viral sequences were enriched and bacterial and archaeal sequences were depleted in viromes, relative to total metagenomes, we performed a series of analyses to compare these two approaches. After quality filtering, total metagenomes yielded an average of 8,741,015 paired reads per library for April samples and 14,551,631 paired reads for August samples, while viromes yielded an average of 9,519,518 and 5,770,419 paired reads in April and August, respectively (Fig. 1A and Supplementary Table 2). Viromes displayed a significant depletion of bacterial and archaeal sequences, as evidenced by fewer reads classified as 16S rRNA gene fragments: 0.006% of virome reads, compared to 0.042% of reads in total metagenomes (Fig. 1B). Moreover, taxonomic classification of the recovered 16S rRNA gene reads revealed clear differences in the microbial profiles associated with each approach: total metagenomes were significantly enriched in Acidobacteria, Actinobacteria, Firmicutes, and Thaumarchaeota, whereas viromes were significantly enriched in Armatimonadetes, Saccharibacteria, and Parcubacteria (Supplementary Fig. 2A, B). These last two taxa belong to the candidate phyla radiation and are typified by small cells [59,60,61], which would be more likely to pass through the 0.22-um filter that we used for viral particle purification [37, 62]. Although we acknowledge that taxon-specific differences in 16S rRNA gene copy numbers could theoretically account for some of the observed differences in absolute numbers of reads assigned to 16S rRNA genes between viromes and total metagenomes [63], in the context of subsequent analyses (see below), the most parsimonious interpretation is that both the abundances and types of bacterial and archaeal genomic content differed between the two datasets.

To assess differences in sequence complexity between the two profiling methods, we calculated the k-mer frequency spectrum for each library (Fig. 1C). Relative to viromes, total metagenomes displayed an increased number of singletons (k-mers observed only once) and an overall tendency toward lower k-mer occurrences, indicating that size-fractionating our soil communities reduced sequence complexity. These differences in sequence complexity translated into notable contrasts in the quality of de novo assemblies obtained from individual libraries (Fig. 1D), while viromes yielded 800 Mbp of assembled sequences across 169,421 contigs (250 Mbp assembled in ≥10 Kbp contigs), total metagenomes produced only 65 Mbp across 22,951 contigs (1.5 Mbp assembled in ≥10 Kbp contigs). The improved assembly quality from the viromes was despite lower sequencing throughput relative to total metagenomes, particularly for the August samples (Fig. 1A). Using DeepVirFinder [48] and VirSorter [47] to mine assemblies for viral contigs, we found that 52.4% of virome contigs and only 2.2% of total metagenome contigs were identified as viral. Together, these results show that our laboratory methods for removing contamination from cells and free DNA reduced genomic signatures from cellular organisms, substantially improved sequence assembly, and successfully enriched the viral signal in soil viromes relative to total metagenomes.

Viromes facilitate exploration of the rare virosphere

To remove redundancy in our assemblies, we clustered all 192,372 contigs into a set of 105,909 representative contigs (global identity threshold = 0.95). Following current standards to define viral populations (vOTUs) [51, 52], we then screened all nonredundant ≥10 Kbp contigs for viral signatures. We identified 4065 vOTUs with a median sequence length of 17,870 bp (max = 259,025 bp) and a median gene content of 27 predicted ORFs (max = 421 ORFs). To profile the viral communities in our samples, we mapped reads against this database of nonredundant vOTU sequences (≥90% average nucleotide identity, ≥75% coverage over the length of the contig). On average, 0.04% of total metagenomic reads and 23.4% of viromic reads were mapped to vOTUs (Supplementary Fig. 3A). One August virome sample (CS-H) had particularly low sequencing throughput and low vOTU recovery (Fig. 1A and Supplementary Fig. 3B) and was discarded from downstream analyses.

In total, 2961 vOTUs were detected through read mapping in at least one sample. Of these, 2864 were exclusively found in viromes, 94 in both viromes and total metagenomes, and three in total metagenomes alone. Thus, viromes were able to recover 30 times as many viral populations as total metagenomes, even when vOTUs assembled from viromes were part of the reference set for read mapping. Notably, the three vOTUs exclusively detected in total metagenomes were only present in one metagenome from April that did not have a successful paired virome (Supplementary Fig. 4). Considering that all other vOTUs detected in total metagenomes were detected in at least one virome, it seems possible that the corresponding virome could have contained these vOTUs if sequencing had been successful. Consistent with capturing a representative amount of viral diversity from the viromes but not total metagenomes, our sampling effort was sufficient to approach a richness asymptote in vOTU accumulation curves derived from viromes but not total metagenomes (Fig. 2A).

To examine the distribution of vOTUs along the abundance-occupancy spectrum, we compared mean relative abundances of vOTUs against the number of samples in which each vOTU was detected. Given the contrasting experimental conditions between the April and August collections, we performed this analysis within each time point. In viromes, highly abundant vOTUs tended to be recovered in the majority of samples (i.e., they displayed high occupancies), while rare vOTUs were typically recovered in only a few samples (Fig. 2B, C and Supplementary Fig. 5A), a trend usually observed in microbial communities [64]. Furthermore, more than 30% of vOTUs were found in all sampled plots, indicating the presence of a sizable core virosphere distributed throughout the field. In contrast, the distribution of 16S rRNA gene OTUs in viromes leaned toward lower occupancies (Supplementary Fig. 5B) as expected from the significant depletion of cellular genomes upon size fractionation (Fig. 1B). On the other hand, more than 80% of vOTUs in total metagenomes were detected only once (Supplementary Fig. 5C), despite the widespread distribution displayed by the 16S rRNA gene OTUs identified in the same samples (Supplementary Fig. 5D), suggesting a sparse recovery of viral diversity compared to a more complete recovery of bacterial and archaeal diversity in total metagenomes.

Inspecting the abundance-occupancy patterns for the 94 vOTUs detected in both viromes and total metagenomes revealed that vOTUs recovered from total metagenomes were among the most abundant and ubiquitous in virome profiles (Fig. 2B), indicating that total soil metagenomes were more likely to miss the rare virosphere. Notably, comparing the relative abundances of vOTUs across paired total metagenomes and viromes showed that their abundance-based ranks were not always preserved (Supplementary Fig. 4). While this discrepancy could stem from methodological challenges associated with virome preparations (e.g., differential adsorption of viruses to the soil matrix could have impacted their resuspension and recovery, therefore affecting the relative abundances of the associated vOTUs), it is also likely that total metagenomes were more susceptible to subsampling biases as evidenced by the sparse and inconsistent recovery of vOTUs exhibited by this profiling method (Supplementary Fig. 5C).

Viromes reveal a diverse taxonomic landscape

To examine the taxonomic spread covered by our vOTUs, we compared them against the RefSeq prokaryotic virus database using vConTACT2, a network-based method to classify viral contigs [49]. Under this approach, vOTUs are grouped by shared predicted protein content into taxonomically informative viral clusters (VCs) that approximate viral genera. Of the 2961 vOTUs, 1712 were confidently assigned to VCs, while the rest were only weakly connected to other clusters (outliers, 784 vOTUs) or shared no genus-level predicted protein content with any other contigs (singletons, 465 vOTUs) (Supplementary Table 3). Only 130 vOTUs were grouped with RefSeq genomes, indicating that this dataset has substantially expanded known viral taxonomic diversity (Fig. 3A). Subsetting the vOTUs detected by each profiling method showed that viromes captured a more taxonomically diverse set of viruses: 1711 vOTUs detected in viromes were assigned to 533 VCs, while 68 vOTUs detected in total metagenomes (67 of which were also detected in viromes) were assigned to 54 VCs (53 of which were also detected in viromes). Thus, any potential biases in the types of viruses recovered through soil viromics (e.g., through preferential recovery of certain viral taxonomic groups from the soil matrix) were not immediately obvious. Any such biases were either eclipsed by the much greater taxonomic diversity of viruses recovered in viromes, relative to total metagenomes, and/or they also apply to total metagenomes, at least for the viruses and soils examined here.

Most of the 130 vOTUs clustered with RefSeq viral genomes could be taxonomically classified at the family level (Fig. 3B). Podoviridae was the most highly represented family, followed by Siphoviridae and Myoviridae. Myoviridae were only detected in viromes, not total metagenomes, further confirming that viromes do not seem to exclude viral groups relative to total metagenomes—if anything, the opposite seems to be true. Among the Siphoviridae clusters, we could further identify three vOTUs as belonging to the genus Decurrovirus, which are phages of Arthrobacter, a genus of Actinobacteria common in soil [65,66,67]. Because the genome network was highly structured by host taxonomy (Fig. 3A, [68]), we used consistent host signatures among RefSeq viruses in the same VC to assign putative hosts to vOTUs in VCs shared with RefSeq genomes. Most such vOTUs were putatively assigned to Proteobacteria, Actinobacteria, or Bacteroidetes hosts, and a few were linked to Firmicutes. Interestingly, these bacterial phyla were among the most abundant taxa in the 16S rRNA gene profiles derived from the total metagenomes from these soils (Supplementary Fig. 2A).

Although soil viromes and total metagenomes have been compared [62] and their presumed advantages and disadvantages have been reviewed [10], here a comprehensive comparison of results from both profiling approaches applied to the same samples showed that soil viromes recover richer (Fig. 2A) and more taxonomically diverse (Fig. 3) soil viral communities than total metagenomes.

Compositional patterns of agricultural soil viral communities and their ecological drivers

Since viromes vastly outperformed total metagenomes in capturing the viral diversity in our samples, we focused on viromes to explore the compositional relationships among viral communities. To assess the impact of each individual experimental factor on beta diversity, we performed separate permutational multivariate analyses of variance (PERMANOVA) on Bray–Curtis dissimilarities (Supplementary Table 4). Collection time point had a significant effect (R² = 0.50, p = 0.001), but biochar (R² = 0.19, p = 0.58) and nitrogen (R² = 0.12, p = 0.75) treatments did not (only samples from the August time points, after nitrogen amendments, were considered for the nitrogen analysis). Additionally, to determine if the location of each sampled plot had an impact on community composition, we tested the effect of plot position along the West–East (W–E) and South–North (S–N) axes of the field (Supplementary Table 4). Viral communities displayed a significant spatial gradient along the W–E axis (R² = 0.17, p = 0.046) but not the S–N axis (R² = 0.10, p = 0.20). Given the significant spatiotemporal structuring in our samples, we performed an additional PERMANOVA to examine the effect of biochar while accounting for these factors (Supplementary Table 5) and detected a significant effect (R² = 0.19, p = 0.012) only when both collection time point and W–E gradient were part of the model. We did not detect a significant effect of nitrogen treatment, even after accounting for the W–E gradient in the August samples (Supplementary Table 5).

To assess whether the bacterial and archaeal communities displayed similar compositional trends and could therefore potentially explain patterns in viral community composition, we attempted to generate metagenome assembled genomes (MAGs) from our total metagenomes. However, the low quality of total metagenomic assemblies (Fig. 1D) precluded MAG reconstruction (19 MAGs with a median completeness of 30.3), so instead we used 16S rRNA gene profiles recovered from total metagenomes (Supplementary Fig. 2A). Although 16S rRNA genes accounted for <0.05% of the reads in our total metagenomes (Fig. 1B), 573 bacterial and archaeal 16S rRNA gene OTUs were detected, and accumulation curves approached richness asymptotes (Supplementary Fig. 2C), suggesting that enough microbial diversity was recovered to justify further analyses. PERMANOVA on Bray–Curtis dissimilarities revealed that only collection time point significantly correlated with bacterial and archaeal community composition, while spatial location, biochar treatment, and nitrogen fertilizer concentration did not (Supplementary Tables 4 and 5).

Viral and microbial communities display coupled temporal dynamics

Compositional differences between April and August samples were significant for both viral and microbial communities (Fig. 4A), and a Mantel test revealed that the Bray–Curtis dissimilarity matrices of viral and microbial communities were significantly correlated (r = 0.59, p = 0.0003). Unsurprisingly, due to the reliance of viruses on their hosts for replication, viral and microbial communities have been previously observed to correlate in soil [14, 18], and results here suggest that viruses and their bacterial and archaeal hosts had coupled temporal responses to the same variables and/or to each other.

To identify individual vOTUs that displayed temporal dynamics, we performed a differential abundance analysis between collection time points. Of 2958 vOTUs detected in viromes, 1062 were enriched in April samples, and 681 were enriched in August samples (Supplementary Table 6). The summed relative abundances of these vOTUs accounted for up to 63.3% of the viral communities (Fig. 4B). No clear taxonomic trend could be detected, as very few (92 of 1743) differentially abundant vOTUs could be taxonomically classified, and all identified viral families and putative host phyla were proportionally represented across April- and August-enriched vOTUs (Supplementary Fig. 6A, B).

To identify microbial taxa associated with compositional shifts between time points, we performed a similar differential abundance analysis on the 16S rRNA gene OTU profiles. Of 573 taxa, 38 were enriched in April, and 20 were enriched in August (Supplementary Table 7). This relatively small number of 16S rRNA gene OTUs accounted for 66.3% of the total microbial community abundance (Fig. 4B). Taxa enriched in April encompassed a diverse set of microorganisms, with members of Bacteroidetes, Alphaproteobacteria, Gammaproteobacteria, and Acidobacteria being the most prevalent, while taxa enriched in August were exclusively members of the Actinobacteria or Alphaproteobacteria (Supplementary Fig. 6C). Interestingly, a recent study on greenhouse-grown tomatoes [69] showed similar taxonomic enrichment and depletion patterns in rhizosphere communities relative to bulk soil communities, including enrichment for Actinobacteria and Alphaproteobacteria in the rhizosphere, suggesting that root microbiome recruitment could explain some of the compositional changes between preplanting and harvest time points. These results are consistent with many previous studies that have demonstrated rhizosphere impacts on the diversity and composition of microbial communities [70,71,72,73] and, more recently, on RNA viral communities in Mediterranean grasslands [16].

We also inspected soil properties and nutrient profiles to determine whether they followed the same temporal patterns as viral and microbial communities. Missing colon. Indeed, they did: samples grouped into two distinct clusters based on overall chemical composition, one with all April samples and one with all August samples (Fig. 4C). Testing the impact of collection time point on individual chemical measurements revealed that, relative to April soils, August soils exhibited a significant increase in ammonium, nitrate, total nitrogen, sodium, and sulfate concentrations and a significant decrease in soil moisture and calcium content (ANOVA, Supplementary Table 8). Nitrogen amendments, which were applied between sampling time points, have been shown to shift the composition of microbial communities associated with tomato rhizospheres toward an Actinobacteria-enriched state [23], suggesting that the increased abundance of this phylum in August samples may be related to an increase in nitrogen availability, in addition to the presence of tomato roots. Furthermore, a recent survey of rice paddies revealed that nitrogen amendments influenced the absolute abundances of viruses and bacteria [74], indicating that a coupled response by viral and microbial communities could have been triggered by agricultural nitrogen inputs. Considering that soil moisture availability is one of the main factors influencing the composition of soil bacterial communities [75] and that soil desiccation has been linked to viral inactivation [24], it is likely that soil moisture also contributed to viral and bacterial community differences between time points. Similar correlations between viral community composition and soil moisture content have been observed in thawing permafrost peatlands [14].

Altogether, these results indicate that viral and microbial communities displayed strong temporal dynamics that were consistent with responses to changes in both biotic (plant) and abiotic characteristics of soils. Still, we acknowledge that the magnitude of this temporal shift could have been amplified by biases associated with library preparation: while April samples were constructed using a standard ligation-based approach, August samples were inadvertently constructed with a tagmentation protocol due to changes in methodology at the sequencing facility (see “Methods”). Given that nonrandom transposition might skew the compositions of viromic [76] and metagenomic [77] profiles, it is possible that some of the detected temporal differences in our study could stem from technical artifacts. However, considering the consistency of some of our results with previous studies (e.g. specific taxonomic response patterns associated with rhizosphere assembly and nitrogen fertilization), it is likely that, overall, the observed compositional differences represent true ecological patterns over time.

Viral communities but not microbial communities were spatially structured in the West–East direction across an agricultural field

To visualize the effect of plot position on the overall structure of viral communities, we performed a principal coordinate analysis on Bray–Curtis dissimilarities calculated across the virome-derived vOTU profiles. While the first axis captured the differences between collection time points, the second axis revealed a compositional transition from the westernmost (W) to the easternmost (E) end of the field (Fig. 5A, B). As expected from the PERMANOVA results (Supplementary Tables 4 and 5), the unconstrained ordination did not display any discernible patterns based on biochar treatment or nitrogen fertilization regime (Supplementary Fig. 7A, B). Comparing the pairwise Bray–Curtis dissimilarities in viral community composition against the corresponding pairwise W–E field distances further confirmed a significant correlation between viral community composition and W–E position (Fig. 5C). This pattern was not recapitulated along the South–North direction. We found 1035 vOTUs displaying a W–E gradient in their relative abundances across the field: 460 with increasing abundances from east to west and 575 with increasing abundances from west to east (Fig. 5D and Supplementary Table 9). Furthermore, examining the overlap with vOTUs affected by collection time point (Fig. 4B) revealed that more than half of the W–E spatially structured vOTUs were temporally dynamic (Supplementary Table 9 and Fig. 5D), suggesting that this W–E structuring of viral communities occurred on time scales captured by this study. In other words, the vOTUs patterned by a W–E gradient in April were generally distinct from those exhibiting similar W–E abundance trends in August, presumably indicating either repeated responses to the same underlying variables and/or repeated patterns of dispersal over time.

We next explored whether this W–E gradient was observed in the bacterial and archaeal communities. As expected from the PERMANOVA results (Supplementary Tables 4 and 5), there was no significant correlation between microbial community composition and spatial distance across the W–E field axis (Supplementary Fig. 8A, B). Thus, while the overarching temporal differences in viral and microbial communities were related, the W–E spatial distribution of viral and microbial communities would seem to have resulted from at least partially decoupled assembly processes. This result is consistent with Mantel tests between viral and microbial community composition, which were significant overall (i.e., driven by temporal separation, see above), but were not significant within each collection time point (April: r = 0.12, p = 0.27; August: r = 0.27, p = 0.15). We also analyzed the soil nutrient profiles and did not detect a significant effect of W–E position on any of the measured chemical or environmental properties (Supplementary Table 8).

Given the strong spatiotemporal dynamics displayed by viral communities, we examined the effect of biochar treatment on overall viral community composition through a partial canonical analysis of principal coordinates (CAP) that removed the variance due to collection time point and W–E position. This approach revealed that biochar-amended soils clustered separately from untreated soils and that viral communities associated with different types of biochar were compositionally distinct, a pattern driven by 43 differentially abundant vOTUs (Supplementary Fig. 9A, B and Supplementary Table 10). Similar trends were observed in the 16S rRNA gene profiles of total metagenomes despite the effect of biochar on overall microbial community composition being not significant (Supplementary Fig. 9C and Supplementary Tables 4 and 5). While previous studies have detected clear compositional shifts in the bacterial communities of agricultural fields upon biochar additions [78,79,80], no consistent pattern has emerged among responsive taxa across different biochar amendments and soil types [81]. Furthermore, while biochar was spread evenly throughout the field in these previous studies, the subsurface banding approach followed in this study constrained the biochar to a narrow area along each bed (Supplementary Fig. 1A). This could explain why none of the measured chemical soil properties was significantly influenced by biochar treatment (Supplementary Table 8), as soil samples were collected 15 cm away from the biochar application zone (Supplementary Fig. 1A). Overall, these results suggest that biochar amendments may have impacted both viral and microbial community composition but that temporal differences (and, for viral communities, W–E spatial differences) were far more pronounced. Thus, future studies that seek to investigate the impacts of soil amendments on viral communities may benefit from rigorous spatiotemporal replication in the study design.

While spatial differentiation has been previously reported for soil viral communities, prior studies considered larger areas with contrasting soil properties. A RAPD DNA study of soil viral diversity along a 4 km land-use transect (including forest, pasture, and cropland systems) found significant differences in viral community composition along the transect [82]. Similarly, a metagenomic characterization of peatland soils found distinct viral communities in three areas comprising different habitats along a permafrost gradient within a ~12,500 m² area [14]. Although the spatial patterns in viral community composition detected in these studies correlated with differences in various environmental parameters, none of the soil properties that we measured in this study were able to explain the observed W–E spatial gradient. Moreover, we found no evidence of W–E spatial structuring in the coexisting bacterial and archaeal communities, suggesting that viruses were uniquely affected by an unknown underlying factor.

One possible explanation for the observed W–E spatial structuring of viral communities could be unidentified legacy effects from previous growing seasons that may have created an unrecognized gradient across the agricultural field (though we note that the field was tilled in both S–N and W–E directions shortly before biochar addition, likely dampening any potential legacy effects, and the field was fallow for at least one year prior to the start of our study). It has been proposed that rates of viral decay might be much slower than viral production in some soils [83], such that an accumulation of older viruses might explain differences in the W–E spatial structuring of viral, relative to microbial communities. However, if this were the case, we would expect most of the W–E spatially structured vOTUs to persist over time and exhibit the same abundance patterns in the April and August samples, yet more than half of these vOTUs were also differentially abundant between time points. Another potential factor that could have contributed to this W–E spatial differentiation is dispersal, consistent with the location of the sampled field next to a gravel road (Supplementary Fig. 1A). We speculate that vehicular traffic could have transported dust from the road to the field and/or provided vibrations to induce shifts in soil pore structure and hydraulic connectivity, both of which would presumably have been more pronounced closer to the source (the road). Additionally, contrasting conditions at the W–E borders of the field (a road on one side and more tomato fields on the other, Supplementary Fig. 1A) could also have led to edge effects. While little is known about the impact of edge effects on soil viruses, other biotic and abiotic properties of agricultural soils have been shown to be influenced by their proximity to the field border [84, 85]. Future studies could be designed to better disentangle the complex interactions of production, decay, selection, and dispersal on soil viral and microbial community composition.

Source: Ecology - nature.com