key: cord-0273427-hm49cfkm
authors: Santos-Medellín, Christian; Estera-Molina, Katerina; Yuan, Mengting; Pett-Ridge, Jennifer; Firestone, Mary K.; Emerson, Joanne B.
title: Spatial turnover of soil viral populations and genotypes overlain by cohesive responses to moisture in grasslands
date: 2022-03-24
journal: bioRxiv
DOI: 10.1101/2022.03.24.485562
sha: d79ab9f36abed94381ea9394017028ce5618c1ef
doc_id: 273427
cord_uid: hm49cfkm

Although soil viral abundance, diversity, and potential roles in microbial community dynamics and biogeochemical cycling are beginning to be appreciated1–5, little is known about the patterns and drivers of soil viral community composition that underlie their contributions to terrestrial ecology. Here, we analyzed 43 soil viromes from a precipitation manipulation experiment in a Mediterranean grassland in California, USA. We recovered 5,315 viral population sequences (vOTUs), and viral community composition exhibited a highly significant distance-decay relationship within the 18 m long field. This pattern was recapitulated in the microheterogeneity of 130 prevalent vOTUs (detected in >=90% of the viromes), which tended to exhibit significant negative correlations between genomic similarity of their predominant allelic variants and distance. Although spatial turnover was also observed in the bacterial and archaeal communities from the same soils, the signal was dampened relative to the viromes, suggesting differences in assembly drivers at local scales for viruses and their microbial hosts and/or differences in the temporal scales captured by viromes and total DNA. Despite the overwhelming spatial signal, vOTUs responsive to a decrease in soil moisture were significantly enriched in a predicted protein-sharing subnetwork of 326 vOTUs linked to 191 known actinobacteriophages, suggesting a genomically cohesive viral response to soil moisture evocative of environmental filtering, potentially by way of actinobacterial hosts. Overall, soil viral ecological processes appear to be highly constrained in space and tightly coupled to the heterogeneous, dynamic soil environment and thus fundamentally different from those of their well-mixed and more thoroughly studied marine counterparts.

about the patterns and drivers of soil viral community composition that underlie their 16 contributions to terrestrial ecology. Here, we analyzed 43 soil viromes from a precipitation 17 manipulation experiment in a Mediterranean grassland in California, USA. We recovered 18 5,315 viral population sequences (vOTUs), and viral community composition exhibited a 19 highly significant distance-decay relationship within the 18 m long field. This pattern was 20 recapitulated in the microheterogeneity of 130 prevalent vOTUs (detected in >=90% of 21 the viromes), which tended to exhibit significant negative correlations between genomic 22 similarity of their predominant allelic variants and distance. Although spatial turnover was 23 2 also observed in the bacterial and archaeal communities from the same soils, the signal 24 was dampened relative to the viromes, suggesting differences in assembly drivers at local 25 scales for viruses and their microbial hosts and/or differences in the temporal scales 26 captured by viromes and total DNA. Despite the overwhelming spatial signal, vOTUs 27 responsive to a decrease in soil moisture were significantly enriched in a predicted 28 protein-sharing subnetwork of 326 vOTUs linked to 191 known actinobacteriophages, 29 suggesting a genomically cohesive viral response to soil moisture evocative of 30 environmental filtering, potentially by way of actinobacterial hosts. Overall, soil viral 31 ecological processes appear to be highly constrained in space and tightly coupled to the 32 heterogeneous, dynamic soil environment and thus fundamentally different from those of 33 their well-mixed and more thoroughly studied marine counterparts.

Main text 36 With an estimated area of 52.5 million square kilometers 6 , grasslands are major 37 contributors to the cycling 7 and storage 8 of soil organic carbon at a global scale. Soil 38 microorganisms play key roles in these biogeochemical processes 9,10 , and, by infecting 39 soil microbiota 11 , viruses likely have substantial direct and indirect impacts on the 40 resulting carbon dynamics 12 . More generally, the potential importance of viruses in 41 soils 1, 2, 13, 14 , together with their measured high abundance (10 7 to 10 10 virus-like particles 42 per gram of soil 1 ) and improvements in our ability to sequence and track soil viral 43 genomes 12, 15 , has led to a renewed flurry of investigations into soil viral diversity and 44 ecology 3-5,16-21 . Yet, despite a new appreciation for the vast diversity of soil viruses 3-5,16-45 18 , little is known about the factors that govern soil viral community assembly. (Figure 1b) . A significant negative correlation 67 between Bray-Curtis similarity and spatial distance between plots (Figure 1c ) further 68 revealed that distance-decay relationships were a key driver of viral community 69 4 composition. These trends were reinforced by substantial differences in individual vOTU 70 detection patterns: of 5,135 vOTUs, 50% were detected in 9 or fewer of the 43 samples 71 (Supplementary Figure 2b) . Moreover, the percentage of vOTUs shared between pairs 72 of viromes declined steeply as spatial separation increased (Supplementary Figure 2c) . 73 To assess whether the bacterial and archaeal communities displayed similar 74 spatial patterns, we performed 16S rRNA gene amplicon profiling on total DNA extracted 75 from the same soil samples used to generate the viromes. While some spatial structuring 76 of the bacterial and archaeal communities was observed, this pattern was only evident 77 along the fourth and fifth axes in a principal coordinates analysis (Figure 2d) . Similarly, 78 even though spatial distance was significantly negatively correlated with microbial 79 community Bray-Curtis similarity (Figure 2e) , this association was not as pronounced as 80 for the viral communities. In particular, the turnover rate of community similarity over 81 distance (the slope of the distance-decay relationship) for viruses was 5.7 times higher 82 than for bacteria and archaea. These differences in the strength of the spatial patterns 83 between viruses and prokaryotes could be related to differences in the integrated In addition to environmental selection and dispersal, diversification (i.e., the 115 generation of novel genetic variation) can contribute to diversity patterns in microbial 116 communities 30,32,33 . To explore the role of spatial structuring on viral genotypic 117 heterogeneity across the field site, we profiled within-population genomic variation. 118 Briefly, using inStrain 34 , we scanned all mapped reads assigned to individual vOTUs and 119 identified polymorphic sites. Then, we reconstructed sample-specific consensus vOTU 

Viruses in 467 Soil Ecosystems: An Unknown Quantity Within an Unexplored Territory

Viruses in soil: Nano-scale undead drivers of 470 microbial life, biogeochemical turnover and ecosystem functions

Host-linked soil viral ecology along a permafrost thaw gradient

Active virus-host interactions at sub-freezing temperatures in Arctic 475 peat soil

Hidden diversity of soil giant viruses

Food and Agriculture Organization of the United Nations

Biogeochemical cycling in grasslands 480 under climate change

A trade-off between plant and soil carbon storage under elevated 482 CO2

Microbial contributions to climate change 484 through carbon cycle feedbacks

The role of soil microbes in the global 486 carbon cycle: tracking the below-ground microbial processing of plant-derived carbon 487 for manipulating carbon dynamics in agricultural systems

Stable-Isotope-Informed, Genome-Resolved Metagenomics 490 Uncovers Potential Cross-Kingdom Interactions in Rhizosphere Soil

Soil Viruses: A New Hope. mSystems 4

The 'Neglected' Soil Virome -Potential Role and 494

Embracing the unknown: disentangling the complexities of the soil 496 microbiome

Coming-of-Age Characterization of 498

Soil Viruses: A User's Guide to Virus Isolation, Detection within Metagenomes, and 499 Viromics

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

Minnesota peat viromes reveal terrestrial and aquatic niche 504 partitioning for local and global viral populations

Metatranscriptomic reconstruction reveals RNA viruses with the potential to shape 507 carbon cycling in soil

Soil Viruses Are Underexplored Players in Ecosystem Carbon 509 Processing

de-replication

Ultrafast and memory-efficient 606 alignment of short DNA sequences to the human genome

The Sequence Alignment/Map format and SAMtools

Prodigal: prokaryotic gene recognition and translation initiation site 611 identification

Fast and sensitive protein alignment using 613 DIAMOND

Graph clustering by flow simulation

Extension to 'ggplot2'

SortMeRNA: fast and accurate filtering of 617

ribosomal RNAs in metatranscriptomic data

The SILVA ribosomal RNA gene database project: improved data 619 processing and web-based tools

Naive Bayesian classifier for 621 rapid assignment of rRNA sequences into the new bacterial taxonomy

Ribosomal Database Project: data and tools for high throughput 624 rRNA analysis

627 PANDAseq: paired-end assembler for illumina sequences

Search and clustering orders of magnitude faster than BLAST

QIIME allows analysis of high-throughput community 632 sequencing data

APE: Analyses of Phylogenetics and Evolution 636 in R language

Associations between species and groups of sites: 638 indices and statistical inference

Multcomp: simultaneous inference in general parametric models. R 640 package version

Moderated estimation of fold change and 642 dispersion for RNA-seq data with DESeq2

Elegant Graphics for Data Analysis

The igraph software package for complex network 646 research. InterJournal, complex systems 1695

) and lower (red) blocks. Square 653 outlines indicate the rainfall manipulation regime assigned to each plot. Differences in 654 font color are for legibility only. (b,d) Unconstrained analysis of principal coordinates 655 performed on (b) vOTU and (d) 16S rRNA gene OTU Bray-Curtis dissimilarities. Panel 656 (b) displays the first and second axes while panel (d) displays the fourth and fifth axes, 657 as they best captured the spatial structuring in (b) viral and (d) bacterial and archaeal 658 communities. Color reflects the plot from which the sample was collected and 659 corresponds to the gradient palette in panel (a)

Pearson's correlation coefficient (r), the linear regression slope, and the associated P

(measured as the percentage of 673 polymorphic sites in a vOTU sequence) and (b) mean relative abundance within prevalent 674 (≥ 90% occupancy) and non-prevalent (< 90% occupancy) vOTUs (c) Distributions of 675 average nucleotide identities (ANI) for each prevalent vOTU, calculated between pairs of 676 sample-specific vOTU consensus sequences. Each box plot corresponds to a single 677 vOTU, and the y-axis is in rank order (ascending from top to bottom) of the median ANI 678 value for each vOTU. Boxes display the median and interquartile range (IQR), and data 679 points farther than 1.5x IQR from box hinges are plotted as outliers. The heatmap on the 680 right shows the Pearson's correlation coefficients between consensus ANI and spatial 681 distance. Bold outlines indicate a significant P-value (< 0.05) for the correlation after 682 multiple comparisons correction (Holm algorithm). Filled black squares correspond to 683 vOTUs with no variation across samples (i.e., all ANIs were equal to 1). (d) The top 5 684 vOTUs with the most significant correlations (lowest P-values) between consensus ANI 685 and spatial distance

Note that vOTUs are defined in part by sharing >= 95% ANI (see Methods), so within-689 vOTU ANI values will necessarily be >= 95% ANI

695 and (c), samples are grouped along the x-axis by collection time point (T1 and T2) and 696 precipitation regime (100% and 50%). (a) Distribution of scores along the third axis of a 697 principal coordinates analysis performed on vOTU Bray-Curtis dissimilarities. The y-axis 698 label indicates the percentage of total variance explained. The first two axes of the same 699 analysis are shown in Figure 1b. (b) Gravimetric soil moisture contents. Boxes display 700 the median and interquartile range (IQR), and data points farther than 1.5x IQR from box 701 35 hinges are plotted as outliers

Summed mean relative abundances of the sets of vOTUs detected as indicator species 704 differentiating T2-50 communities from the rest of the viromes. Facets distinguish 705 indicator vOTUs that were relatively enriched or depleted, respectively

Node color shows whether a vOTU was an indicator 708 species enriched or depleted in T2-50 samples or not an indicator species (defined by P-709 values below or above 0.05, respectively, from an indicator value permutation test)

Nodes surrounded by squares 713 correspond to vOTUs with a significant overlap in their predicted protein contents with 714 any of 971 RefSeq phage genomes

All such RefSeq phage genomes with significant links to this 716 subnetwork were from phages isolated on Actinobacteria hosts, indicated by tagging 717 vOTU nodes linked to RefSeq actinophages with the letter "A

Squares mark the locations of 724 individual plots, and the numbers indicate the number of circular subplots from which 725 samples were collected within each plot. Color indicates the rainfall manipulation 726 treatment for each plot. (b) Example of a circular subplot from which samples were 727 collected in this study. Plexiglass dividers segmented subplots into two halves, and one 728 half was destructively harvested at each time point

)) indicate periods when the 732 rainfall-excluding shelters were deployed for each treatment. (d) Differential watering 733 events during the months preceding sample collection. Each bar indicates the amount of 734 water added to individual plots through irrigation, based on their assigned watering regime 735 (colored as in panel (a)). (e) Growth patterns of Avena barbata during the 2020 growing 736 season. Each line displays the average height of A. barbata in a single plot

Supplementary Figure 2. (a) Number of vOTUs detected at each occupancy level

Relationship between the percentage of vOTUs shared across pairs of samples and 743 spatial distance between plots. Each point represents a pair of samples and the spatial 744 distance between them was measured as the length of the line connecting the centers of 745 the corresponding plots

Each facet corresponds to a single variable with the y-axis indicating the 757 absolute values of the differences between pairs of samples. In both (a) and (b), each 758 point represents a pair of samples and the spatial distance between them was measured 759 as the length of the line connecting the centers of the corresponding plots. Trend lines 760 display the least squares linear regression models. Statistics correspond to Pearson's 761 correlation coefficient (r) and associated P-value. Variables with a significant correlation 762 in (b) are highlighted in red. Abbreviations displayed in the facet names correspond to 763 ppm = parts per million; meq/100g = milliequivalents per 100 grams of soil

Node color shows the trait 770 assignment for each vOTU, i.e., whether the vOTU was an indicator species enriched or significant overrepresentation of vOTUs enriched (left facet) or 773 depleted (right facet) in T2-50 samples across the network. Each colored point denotes 774 the center of a significant local neighborhood, representing a total of 10 to 94 vOTUs per 775 point. The color gradient indicates the extent of the significance of trait overrepresentation 776 in that neighborhood, with all shades of blue indicating significance and darker shades 777 showing greater significance

Size (upper panel) and trait composition (lower 779 panel) of local neighborhoods with a significant overrepresentation of vOTUs

In the lower panel, each stacked bar plot shows the fraction of indicator 781 vOTUs within a single neighborhood

Each row represents a single vOTU, and its position along the 789 y-axis is determined by its relative enrichment along the field: vOTUs towards the bottom 790 of the y-axis tended to be more enriched on the North-West (left side of the field in Figure 791 1a), while vOTUs towards the top tended to be more enriched on the South-East (right 792 side of the field). The pink tick marks on the left side highlight the indicator vOTUs that 793 were enriched in T2-50 samples and that were part of the subnetwork