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

Viruses impact nearly all organisms on Earth, with ripples of influence in agriculture, health, and biogeochemical processes. However, very little is known about RNA viruses in an environmental context, and even less is known about their diversity and ecology in soil, 1 of the most complex microbial systems. Here, we assembled 48 individual metatranscriptomes from 4 habitats within a planted soil sampled over a 22-d time series: Rhizosphere alone, detritosphere alone, rhizosphere with added root detritus, and unamended soil (4 time points and 3 biological replicates). We resolved the RNA viral community, uncovering a high diversity of viral sequences. We also investigated possible host organisms by analyzing metatranscriptome marker genes. Based on viral phylogeny, much of the diversity was Narnaviridae that may parasitize fungi or Leviviridae, which may infect Proteobacteria. Both host and viral communities appear to be highly dynamic, and rapidly diverged depending on experimental conditions. The viral and host communities were structured based on the presence of root litter. Clear temporal dynamics by Leviviridae and their hosts indicated that viruses were replicating. With this time-resolved analysis, we show that RNA viruses are diverse, abundant, and active in soil. When viral infection causes host cell death, it may mobilize cell carbon in a process that may represent an overlooked component of soil carbon cycling.

We have much to learn about soil, the world’s most diverse and enigmatic microbial habitat. With the advent of metaomics techniques, the diversity, ecology, and impact of bacteria and archaea in soil is being rapidly revealed (1, 2). However, we have only begun to study viruses in soil despite community interest and their probable importance to soil functioning and carbon storage (3). Recent work on double-stranded DNA bacteriophage (phage) in soil has begun to explore their diversity and host interactions (4–6). However, other viral members of the soil community, such as RNA-based phage and eukaryotic RNA viruses, are largely unknown. RNA viruses, those with genomes encoded on RNA, are less studied in environmental contexts because they are not captured in the more common DNA sequencing studies. In some systems, the impact of RNA viruses on community and ecosystem processes has been proposed to rival or exceed the impact of DNA viruses (7, 8). The majority of soil RNA viral work has been single-host focused for agriculturally relevant crops (9) or crop pathogens (10). Environmental RNA virus studies have focused on less complex or more tractable systems, such as marine environments and the human and animal gut (11–14). To our knowledge, no sequencing-based RNA viral community analyses have investigated soil.

Viruses may be major players in biogeochemical cycling. However, much of what is known about viral impacts on elemental cycling comes from aquatic systems. Marine phages can lyse up to one-third of bacteria in ocean waters per day, releasing a huge amount of carbon (15–17). The released components include dissolved organic carbon that is readily metabolized by heterotrophic bacteria, but largely inaccessible to eukaryotic grazers and higher trophic levels (18). This phenomenon, termed “the viral shunt” (18, 19), is thought to sustain up to 55% of heterotrophic bacterial production in marine systems (20). However, some organic particles released through viral lysis aggregate and sink to the deep ocean, where they are sequestered from the atmosphere (21). Most studies investigating viral impacts on carbon cycling have focused on DNA phages, while the extent and contribution of RNA viruses on carbon cycling remains unclear in most ecosystems. Several studies have indicated RNA viruses may outnumber DNA viruses in some instances, a possible hint to their influence in select ecosystems (7, 8, 22). We suspect that, by analogy, lysis of organisms by RNA viruses may represent a large contribution to carbon flow in soils. The process of cell lysis and biomass carbon liberation likely stimulates heterotrophic consumption, with a substantial portion of the carbon lost to the atmosphere as CO2. The liberated cellular debris may be protected from microbial access, through interaction with mineral surfaces or occlusion within soil aggregates (generation of mineral–organic associations) (23, 24). These processes could conceivably result in carbon stabilization and ultimately long-term soil carbon persistence, although in truth we have much to learn about the processes influencing soil carbon stabilization.

Understanding the diversity and ecology of soil viruses may contribute to advances in biotechnology. Viral genomes have been mined for biopesticides and self-assembling nanomaterials (25, 26). Viruses have been proposed, and used, as biocontrol agents for culling invasive organisms, including fire ants and moths (27, 28). Viruses are also being investigated as biocontrol agents for devastating plant pathogens, such as Fusarium sp., Botrytis cinereal, and Rosellinia necatrix (10, 29–34). Novel, environmentally derived viruses may be a source for new biotechnology tools and biocontrol agents.

Here, we used assembled metatranscriptomic data from a California annual grassland soil to reconstruct RNA viral genome sequences and constrain their possible hosts. We searched the assembled sequences for the RNA-dependent RNA polymerase (RdRp), a conserved gene found in all RNA viruses that lack a DNA stage (generally referred to as Baltimore type III–V RNA viruses). Although recent work tracing the deep evolutionary history of RNA viruses has merged viral families into supergroups (35–37), our focus was not on resolving deep phylogenetic viral branches. Rather, our objectives were to investigate fine-scale diversity of RNA viruses, their putative hosts, and their possible functional importance in soil. Thus, we relied heavily on Shi et al. (36), the most comprehensive metatranscriptomic study of RNA virus diversity so far, for taxonomic classification and to delineate viruses from previously undocumented clades.

The diversity of soil eukaryotes that may serve as hosts for RNA viruses remains largely understudied. Generally, RNA viruses are known to infect fungi, plants, animals and the many clades of single-celled eukaryotes, in addition to some Proteobacteria. Soil eukaryotic studies have relied heavily on primer-based sequencing or visual classification, methods that can impart biases and miss novel organisms. Using the genomic information contained in our assembled metatranscriptomes, we identified many clades of eukaryotes without reliance on primers or microscopy. In contrast to many other environmental viral studies, which sequence extracted viral particles, we sequenced whole samples, which included viral genomes and possible host transcripts. This method allowed us to explore both extracellular viruses that would be sequenced in the standard method and viruses within host cells, such as during the infection process or latent infections. We tracked both viral and eukaryotic communities in key soil habitats: The rhizosphere (soil influenced by the root), detritosphere (soil influenced by decaying particulate organic matter; bulk+litter), a combination of the 2 (rhizosphere+litter), and unamended soil (bulk) over time. Using a relatively short sampling timescale (3, 6, 12, and 22 d) allowed us to investigate viral and host community dynamics. Our research places direct constraints on the timescales for virus dynamics and provides a genomic view of RNA viruses in soil.

Results

Experimental Design.

This analysis utilized data generated by a previous study on microbial niche differentiation in soil habitats (38). Briefly, wild oat (Avena fatua), an annual grass common in Mediterranean climates, was grown in microcosms with a sidecar, allowing us to track root age and sample directly from the rhizosphere (39). The experimental set-up used microcosms, half of which included soil mixed with dried ground A. fatua root litter and the other half contained soil without litter amendment. All microcosms contained bulk soil bags, which excluded roots, providing control soil uninfluenced by the root. Once A. fatua was mature, roots were allowed to enter the sidecar and the growth of individual roots was tracked. We destructively harvested rhizosphere soil (and paired bulk soil) that had been in contact with the root for 3, 6, 12, and 22 d. In total, we sampled paired rhizosphere and bulk samples from 4 time points with 2 treatments (with and without litter), with 3 biological replicates, for a total of 48 samples for metatranscriptome sequencing. Using rRNA-depleted RNA, we sequenced a total of 408 Gbp (average of 8.7 Gbp per sample).

Eukaryotic RNA Viruses.

We used profile hidden Markov models (HMM) to search our assembled metatranscriptomes and found a total of 3,884 unique viral RdRp sequences (dereplicated at 99% amino acid sequence identity; AAI). This includes 1,350 predicted RNA bacteriophage (phage) viruses that infect bacteria, and 2,534 predicted viruses that infect eukaryotes.

Our eukaryotic viruses group into 15 major clades that span the majority of known viral diversity (Fig. 1). Many were included into the supergroup-like clades defined by Shi et al. (36). For the remainder, we constructed phylogenetic trees to define additional viral families.

Fig. 1.

Phylogenetic trees representing clades of RNA viruses (based on RdRp) identified in our California annual grassland experimental soil. Within each tree, the RdRp sequences we identified are colored purple and previously described sequences are in pink. Trees are all midpoint rooted. Trees with predicted fungal infecting clades are presented in .

Overall, in trees that include both existing and newly generated sequences, we noted a strong grouping of our RNA viral sequences into “fans of diversity,” much like the seed head of a dandelion. Many of these fans included a single reference sequence previously used to propose a new viral family. For example, in the Hepe-Virga, sequences grouped with the Agaricus bisporus virus 16 [proposed family Ambsetviridae (40)], in the Astro-Poty clade sequences group with Bufivirus UC1 from wastewater (37) and in the Partiti-Picobirna, sequences cluster with Purpureocillium lilacinum nonsegmented virus 1 (41). We substantially expanded the Barnaviridae (Lueto-Sobemo superfamily) and Mymonaviridae (Tombus-noda superfamily) families (42, 43), which replicate in fungi, and the newly proposed Zhaovirus and Qinvirus families, which are found in invertebrates (36).

We predict that fungi are the most common hosts for many of our reported RNA viruses (mycoviruses) (). We are most confident when they fall into the Barnaviridae, Megabirnaviridae, Quadriviridae, and mitoviruses, groups currently thought to only infect fungi (44, 45). The most frequently encountered virus in our dataset, accounting for over 50% of the eukaryotic viral strains identified, came from the mitovirus genus in the Narnaviridae family. Mitoviruses are linear single-stranded RNA viruses that replicate within fungal mitochondria and spread vertically through spores and horizontally through hyphal fusion (33, 46). We suspect that most of the mitoviruses we detect are infecting fungi because, although some mitoviral sequences have been found integrated into the genomes of plants, they are frequently truncated and not transcribed (46, 47). Mitoviruses were also the most abundant viral clade in every sample ().

Until recently, the Narnaviridae group, which includes mitoviruses, were thought to only encode an RdRp, but a recent discovery suggests some narnaviruses (only distantly related to mitoviruses) encode additional proteins, including capsids and helicases (36). We identified some mitovirus genomes with possible additional genes, which could play a part in infection efficiency or infection of new hosts. However, the majority of mitoviruses we identified contained only a single RdRp gene. We predicted several additional proteins on some near-complete mitoviral genomes (). These putative genes are small (average 79 amino acids), and functions could not be predicted for them. Sometimes the additional genes are predicted to be transcribed in the same direction as RdRp and in other cases they are predicted to be transcribed in the opposite direction ().

We reconstructed many sequences for viruses that likely infect eukaryotes other than fungi. For the Picorna–Calici, hosts are likely vertebrates, insects, algae, and plants, based on viral phylogenetic placements (48). In the Tombus-Noda tree, many of the RdRps group with umbraviruses, well-recognized plant viruses. Other RdRp sequences group with sequences from complex environmental samples where the hosts are unknown.

RNA viruses from some previously defined superfamilies were conspicuously absent. We did not identify any soil RNA viruses belonging in the Nidovirales-like, Reoviridae-like or Orthomyxoviridae-like superfamilies. The absence of Reoviridae is interesting, as these RNA viruses can infect fungi (49). We also did not confidently identify any Ortervirales, retroviruses that replicate with a DNA intermediate. We identified many reverse transcriptases, but none of the corresponding scaffolds encoded capsid proteins, so the sequences may be retrotransposons or fragments of retroviruses.

We classified 2 clades within the Bunya-Arena viral superfamily as novel, given that they exhibit sequence divergence comparable to that which separates known RNA viral families (). In both cases, the scaffolds only encode the polymerase with a Bunyavirus RNA-dependent RNA polymerase domain. This genome structure is shared by Hubei myriapoda virus 6 (36) and Ixodes scapularis-associated virus 3 (50).

Eukaryotic Hosts.

As viruses replicate within their host, changes in virus abundance levels implicate introduction of new vectors or hosts, shifts in infected and uninfected host abundance levels, changes in infected host physiology, or changes in host susceptibility to infection. To better understand the diversity and ecology of RNA viral hosts and carriers in soil, we used the mitochondrially encoded cytochrome C oxidase subunit 1 (Cox1) gene as a marker to define the eukaryotic populations present (51–56). However, mitochondrial transcription is determined by metabolic activity of a cell (57) and can be impacted by the switch to a symbiotic state [e.g., in arbuscular mycorrhizal fungi (58)] or nutrient sensing and hormonal signals [e.g., in animals (59)]. Thus, the number of reads mapped to the Cox1 genes is determined by the number of organisms present and the organism’s activity.

We identified 726 eukaryotic Cox1 sequences and clustered them at 98% AAI to approximate a species-level view. The eukaryotes we discovered span the diversity of known soil organisms. Not surprisingly, sequences from A. fatua were the most abundant Cox1 transcripts in many samples. In some samples, the most abundant Cox1 were from an Enchytraeidae-related worm or an Amoebozoa sp. Other samples were dominated by a mixture of fungi, Amoebozoa, Viridiplantae, and unknown eukaryotes.

Amoebozoa were the most diverse eukaryotic clade in this dataset, >25% of the identified eukaryotes, followed by species of fungi and unknown eukaryotes. The many unknown eukaryotes reflect the lack of environmental Cox1 sequences from micro- and mesoeukaryotes in public databases. For this reason, we reconstructed 18S rRNA gene sequences for classification, as they are better represented than Cox1 in public databases. However, the transcripts we analyzed had been depleted in bacterial and plant ribosomal RNA (via Ribo-Zero rRNA Removal kits), which undoubtedly influenced the composition and relative amount of other rRNA sequences. As such, the identified 18S rRNA sequences could not be used for abundance measurement or as a quantitative measure of diversity, only for identification of eukaryotes not identified using COX1. In total we identify 521 distinct species based on 18S rRNA sequences (after clustering at 98% nucleic acid identity to approximate species groups) (Table 1). The results revealed the presence of Centroheliozoa, Malawimonadidae, and Jakobida not identified on the Cox1 analyses and many more species of Heterolobosea, Euglenozoa, and Rhizaria.

Table 1.

Number of eukaryotes identified based on marker genes and significant correlations between eukaryotes and RNA viruses

			Significant cooccurrences
Clade	Cox1	18S rRNA*	Narnaviridae	Bunya-Arena	Tombus-Noda	Luteo-Sobemo	Toti-Chryso	Picorna-Calici	Leviviridae	Parti-Picobirna
Amoebozoa	193	295	14	16	3	2			3	1
Fungi	174	33	15	11	2	1	2	2
Unknown eukaryotes	108	14	1					1
Metazoa	85	28
Stramenopiles	55	17	11	6	3		2	1
Alveolata	54	14
Heterolobosea	11	41
Nucleariidae	11	0	11	2		5				1
Cryptophyta	9	1
Euglenozoa	7	28
Viridiplantae	6	2
Choanoflagellida	5	8	1	3	1
Rhizaria	3	21
Apusozoa	2	1
Alveidia	2	0
Ichthyosporea	1	0
Centroheliozoa	0	11
Malawimonadidae	0	6
Jakobida	0	1

The number of statistically significant correlations is the number of positive cooccurrence edges between a eukaryote and virus derived from the network analysis ().

Note that the absolute values may be skewed by rRNA depletion.

Eukaryotic Virus and Host Ecology.

The presence of root litter shaped the soil’s eukaryotic community structure, as measured by Cox1 gene transcript abundance with A. fatua Cox1 sequences removed (Fig. 2) (PERMANOVA R2 0.17, P < 9.999e-05). Separation of the communities was evident by the time of the first sampling (3 d) and persisted to the end of the 22-d long experiment ().

Fig. 2.

Nonmetric multidimensional scaling ordination of soil eukaryotic communities in response to 4 soil resource treatments, based on coverage of the Cox1 gene (A. fatua sequences removed) (A), and ordination of eukaryotic RNA viral communities with coverage calculated at the scaffold level (B).

We compared the abundances of eukaryotic Cox1 transcripts in rhizosphere+litter, bulk+litter, and rhizosphere samples to the bulk samples in order to identify species statistically enriched in each case (). The results showed that the presence of root litter and rhizosphere enriched for many Amoebozoa and fungi. However, enrichment patterns indicate that litter had a greater selective force on more individual species than the presence of a growing root.

The eukaryotic RNA viral community was influenced by some of the experimental treatments. Root litter had a significant effect on the eukaryotic RNA viral community (PERMANOVA P < 9.999e-05); the presence of growing roots had no detectable impact (PERMANOVA P < 0.07) (Fig. 2). The distinction between added litter vs. no litter samples was evident within 3 d, implying that the viral and eukaryotic communities changed at similar rates. Differences persisted to the end of the 22-d long experiment. We observed the same patterns even when the Narnaviridae were removed from the dataset, indicating that these numerous, unusual host-bound viruses were not causing the separation of the eukaryotic viral community based on the presence of root litter ().

We identified a smaller number of enriched eukaryotic RNA viruses compared the number of enriched eukaryotes (). This may be due to the substantial heterogeneity in viral abundance patterns, even within replicates. The RNA viruses displayed a higher level of microheterogeneity than the eukaryotes. In contrast, transcripts for individual eukaryotes were more ubiquitous across samples. The viral strains most frequently enriched relative to bulk soil were from the mitoviruses. This is unsurprising, as the abundance of the obligately intermitochondrial mitoviruses is tied to the abundance and activity of the host and fungi responded strongly to the specific treatments.

To explore possible correlations between the hosts and RNA viruses we created a cooccurrence network. The network was constructed from positive cooccurrences between the RNA viruses and hosts (). We noted the total number of edges, statistically significant undirected vertices between eukaryotes and RNA virus in the network, between the possible hosts and RNA viruses (Table 1), some eukaryotes and RNA viruses have multiple edges, which explains the high numbers seen in some less diverse clades, such as the Nucleariidae.

RNA Phage, Potential Hosts, and Ecology.

Two families of RNA viruses are known to infect bacteria, the Cystoviridae, which infect Pseudomonas sp. (60, 61), and the Leviviridae, which infect Gammaproteobacteria and Alphaproteobacteria (62–64). We used a marker gene approach to find RdRps for these RNA phage in the assembled transcripts. After dereplication, we identified 12 Cystoviridae RdRp sequences and 1,338 unique Leviviridae RdRp sequences. This is a significant increase in the known diversity of this group, as there are currently just over 200 Leviviridae RdRp sequences in public databases. Some of our sequences group with Allolevivirus and Levivirus, well-studied genera, but importantly, other novel Leviviridae sequences resolved new clades (Fig. 1).

The specific Leviviridae taxa enriched in bulk soil are relatively phylogenetically novel () and have unconventional genetic architectures based on metagenomic reconstruction. The most abundant Leviviridae sequence in 46 of 48 samples groups into a predominantly novel branch. The near-complete 4,668-bp genome of this abundant RNA phage encodes 4 nonoverlapping genes, as opposed to the frequently overlapping and fewer genes found in other Leviviridae (, Top). Like Enterobacteria phage M, it may also encode a +1 frameshift lysis protein inside its RdRp gene (65). The genome of a related Leviviridae is even more complex, with 5 genes and possibly 2 frameshift lysis proteins (, Bottom). However, additional biochemical work will be needed to resolve the function of these predicted proteins.

We used ribosomal protein S3 (rpS3) phylogeny to identify possible hosts for Leviviridae viruses. All previously documented Leviviridae infect Proteobacteria, so we narrowed our analysis to these bacteria, which represent a possible pool of hosts for the Leviviridae, although other clades cannot be ruled out. In our samples, we identified 355 species of Proteobacteria, which differed in abundance between rhizosphere, bulk, rhizosphere+litter, and bulk+litter treatments (PERMANOVA P < 9.999e-05) () after only 3 d of root growth (). The Leviviridae communities also diverged within 3 d and changed over time (PERMANOVA P < 0.02) (Fig. 3). The communities separated based on the presence or absence of root litter (PERMANOVA P < 9.999e-05), whereas the influence of growing roots was undetectable (PERMANOVA P < 0.09) (Fig. 3).

Fig. 3.

Nonmetric multidimensional scaling ordination of the Leviviridae community identified in an annual grassland soil exposed to 4 resource treatments, colored by time (A) and treatment (B).

Discussion

To our knowledge, no prior study has genomically investigated the ecology and diversity of RNA viruses in soils. Using assembled metatranscriptomes, we uncovered a vast diversity of RNA viruses (>3,000 sequences) from a single grassland soil and discovered several possible unique families. Some viruses grouped phylogenetically with viral families previously proposed based on only a single member, indicating soils may hold much of the diversity of these understudied groups. The viruses we reconstructed came from nearly all known groups of RNA viruses. We identified hundreds of eukaryotes and bacteria that may represent RNA viral hosts or vectors that can pass infections among plants and other organisms. We know that the transfer of viral agents can have agriculturally relevant effects (66–68) but there are likely many unexplored impacts on the soil community and host-mediated process, such as nutrient cycling.

Many of the RNA viruses we identified group phylogenetically with fungal viruses (mycoviruses), and fungi appear to be one of the main RNA viral hosts in this soil. The apparent dominance of mycoviruses in our dataset provides a striking contrast to aquatic systems, where Picornavirales dominate (14). The most diverse mycoviruses in our samples are mitoviruses (Narnaviridae) which, despite their simplicity, can have significant impacts on fungal fitness and physiology (33, 69, 70).

Previous studies of environmental RNA virus ecology have primarily addressed marine systems. The RNA viruses of oceans appear to be dominated by Picornavirales, thought to infect diatoms and other single-cell eukaryotes (7, 14, 71). While our samples did contain Picornavirales, they did not appear to be dominant in the soil we sampled. In addition to Picornavirales, some aquatic systems contain plant, human, and livestock RNA viruses that may have been introduced through runoff; we identified no viruses known to infect humans (14, 72). Some aquatic studies have noted viruses related to RNA fungal viruses, however not at the abundances or diversity we have in our soil (7). Again, in contrast to the diversity of RNA phage that we found in soil, there is only a single report of an RNA phage from an aquatic ecosystem (71). Compared to marine systems, soil virology, especially investigations of RNA viruses, is in its infancy.

In our soil, both eukaryotic RNA viral and eukaryotic communities were impacted by the presence of decaying root litter. By the first sampling time point, the eukaryotic and eukaryote-associated RNA viral communities were different in the presence and absence of litter. Saprophytic fungi appeared to respond favorably to the dead litter biomass. In contrast, dissolved organic compounds exuded by roots may be less accessible to nonfungal soil eukaryotes. Accordingly, roots had little impact on the nonplant eukaryotic community composition. Addition of decaying plant biomass promoted the activity of detritivores, which likely immigrated into the litter from surrounding soil, increased transcription levels, grew, and/or germinated from spores. However, the shifts in the eukaryotic RNA viral community we observed are likely due to proliferation following increased activity of their hosts, although they may have also been transported by vectors.

We conducted a network analysis between the RNA viruses and possible hosts to provide new hypotheses about interacting partners. The correlation between virus and host provides possible associations and fodder for future studies. However, overinterpretation of network analysis results, especially in complex systems, should be avoided (73). Within the RNA virus–host cooccurrence network the fungi had the greatest number of cooccurrence links with the Narnaviridae, which provides another piece of evidence that fungi may be hosts for RNA viruses in this soil. Other eukaryotes had cooccurrence links with the Narnaviridae as well, but it is not possible to tell if these viruses are infecting or simply cooccurring with the eukaryote. The Amoebozoa had the largest number of cooccurrence linkages, including 3 with Leviviridae, which are unlikely to infect the Amoebozoa themselves but may be associated with a cooccurring bacteria. We know of only 1 RNA virus that infects Amoebozoa, a Mononegavirales (Mono-Chu) (74). These results suggest a possible link between Amoebozoa and the viral Bunya-Arena clade. The largest module (cluster of linages) contained a phylogenetically diverse array of eukaryotes and viruses including single virus types connected to multiple clades of hosts; this may represent a group of cooccurring entities that respond to similar environmental or biotic conditions ().

RNA phage, which replicate in bacteria, have received relatively little study of their ecological context. Here, we provide a glimpse into their ecology and diversity in soil. Given that over half of all of the identified RdRp sequences in our dataset were Leviviridae, we encourage additional studies (in other soil types, land uses, tillage practices, and cover vegetations) to assess whether soils generally host a large diversity of Leviviridae. As Leviviridae do not have a known lysogenic life stage, most of their population changes likely reflect replication in their hosts. Since Leviviridae communities and their hosts became distinctly different within the first 3 d of our experiment, we conclude that these RNA viruses infected and replicated within days. Infections over this timescale likely have dramatic effects on host communities, and thus soil ecology.

RNA viruses and RNA phage were heterogeneously distributed across samples, including replicates. In contrast, eukaryotes and bacteria appeared to have more even transcript abundances across samples. In combination, these observations suggest that virus and phage abundances are not solely determined by the presence of hosts, but also factors—such as sporadic blooms or variation in viral resistance levels in the host population —that lead to patchy viral distribution patterns in soil. These patterns may be analogous to those documented in marine virus blooms, although the exact causes of these events are still unknown (13, 75). In addition, viruses have different release rates and persistence times in soil, which may be impacted by a variety of factors, such as soil type, viral strain, host physiology, microsite of release, clay content, and extracellular enzymes. This effect would be most pronounced for viruses with an extracellular life stage and less impactful on obligate intracellular viruses like the Narnaviridae.

In general, the magnitude of viral impacts on the soil carbon cycle is underexplored. Little research has been done on phage-induced bacterial lysis in soil and even less on viral-induced death of fungi and other eukaryotes, which can contain an equal or greater biomass compared to bacteria in some soils (76). In addition to lysis, viruses of fungi can have subtle yet profound effects on fungal biology, pathogenesis, mating success, toxin production, symbiotic relationships, and other physiological effects (77–79). These effects of viruses on fungal biology could have impacts for fungal functioning and thus soil performance as a whole. In our samples, the high diversity and abundance of identified RNA viruses, combined with their dynamic population changes, indicates that there was substantial viral replication. We hypothesize that proliferation of these lytic RNA phage and RNA viruses will have substantial impacts on the form, abundance, and distribution of carbon compounds in soil and on the biology and fitness of soil fungi. For example, lysis of host bacterial cells and viral-induced cell death of eukaryotes will cause release of dissolved low molecular weight carbon compounds. These will likely be quickly consumed by nearby bacteria and much of the carbon respired back to the atmosphere, mostly as CO2. However, as happens when marine snow settles into the deep ocean, a portion of the cellular and viral debris and soluble carbon released may be stabilized in soil. For example, bacterial and eukaryotic lipids and polysaccharides could adhere to mineral surfaces or become occluded within soil aggregates.

Conclusions

By reconstructing soil metatranscriptomes drawn from multiple soil habitats and time points within the context of a controlled experiment we greatly increased the known diversity of RNA viruses. Phylogenetic analyses suggest that fungi are the most common hosts for RNA viruses in the studied grassland soil. While the diversity of hosts for RNA phage remain to be definitively identified, they likely also include Proteobacteria. Shifts in eukaryote, RNA phage, and RNA viral abundances over a few-day period reveal that entire soil communities can rapidly respond to altered resource availability. Our experiments indicate that the form of carbon inputs (root-derived low molecular weight C inputs versus macromolecular carbon compound in litter) may impact eukaryotic and bacterial abundance patterns, and that these in turn may be a major determinant of RNA viral and RNA phage dynamics.

Methods

Wild oat (A. fatua) was grown in microcosms with a clear sidecar designed to allow access and visual tracking of the soil and rhizosphere (39, 80). Experimental soil was collected from Hopland Research and Extension Center (Hopland, CA). The soil is a fine-loamy, mixed, active, mesic Typic Haploxeralf, supporting dominant stands of Avena species (81, 82). Microcosms were filled with soil at field bulk density and seedlings were planted and grown for 6 wk in the main chamber before the start of the experiment. Six days before the start of the experiment the divider separating the main chamber and the sidecar was removed and the sidecar was packed with soil. The litter-amended microcosms received 0.4 g of dried A. fatua root litter mixed with 50 g of Hopland soil. Bulk soil was contained in 1-μm mesh bags embedded into the microcosm, which allowed solutes to pass but not roots. Experimental design and soil edaphic characteristics are fully explained in a companion publication, which provides additional details regarding sample collection protocol, RNA extraction and processing, and sequencing (38); what follows is a brief description of the methods.

Sample Collection.

The ages of individual roots were tracked to collect rhizosphere soil that had been influenced by the root for 3, 6, 12, and 22 d. Three replicate microcosms were destructively harvested for paired rhizosphere and bulk soil. Rhizosphere soil was cut out from the rest of the soil along the edge of the root hair zone (<2 mm from the main root). Root sections and adhering soil were placed immediately in ice-cold Lifeguard Soil Preservation Reagent (MoBio) and soil was agitated off the root by vortexing for 2 min on medium speed and pelleted; bulk soil was treated in the same manner. Pelleted soils were frozen on dry ice and stored at −80 °C.

RNA Extraction.

RNA was extracted from 0.5 g of frozen soil using a phenol-chloroform extraction protocol (83). Then Qiagen AllPrep kits to were used to separate DNA from RNA. RNA was treated with TURBO DNase (Thermo Fisher Scientific) following the manufacturer’s protocol to remove residual DNA, and was concentrated using an ethanol precipitation.

Sequencing.

Metatranscriptome libraries were prepared and sequenced at the Joint Genome Institute. Ribosomal RNA was depleted from 1 µg of total RNA using the Ribo-Zero rRNA Removal Kit (Epicentre) for Plants and Bacteria. RNA was sequenced on the Illumina HiSeq sequencing platform utilizing a TruSeq paired-end cluster kit, v3, and Illumina’s cBot instrument to generate a clustered flowcell for sequencing. Sequencing of the flowcell used a TruSeq SBS sequencing kit, v3, following a 2 × 150 indexed run recipe.

Sequence Analysis.

Reads were trimmed using Sickle (https://github.com/najoshi/sickle) and BBtools (https://sourceforge.net/projects/bbmap/) was used to remove Illumina adapters and trace contaminants. The reads from all 48 samples were assembled individually using IDBA-UD with default settings (84). Genes were predicted using Prodigal in the anonymous mode (85).

To find the host marker genes, Cox1 and rpS3, we used an HMM profile from Pfam (86), Cox1 (PF00115), and HMMs for Bacteria and Archaea (https://github.com/AJProbst/rpS3_trckr) and searched using hmmsearch (-E 0.00001) from the HMMER suite (87). The identified proteins were classified using both National Center for Biotechnology Information (NCBI) blast and by making trees. Once classified, we predicted genes again using Prodigal in the single mode with the appropriate translation table. The Cox1 protein sequences were dereplicated and clustered at 98% AAI representing an estimated species-level designation (88, 89). Assembly errors were found and fixed in the scaffolds containing the Cox1 using ra2.py (https://github.com/christophertbrown/fix_assembly_errors) (90). The bacterial and archaeal rps3 genes was clustered at 99% AAI, representing species-level differences (91), using USEARCH (-cluster_fast) (92). The Cox1 and rpS3 trees were generated using references from NCBI, the protein sequences were then aligned using MAFFT v7.402 (93) with the E-INS-i option on Cipres (94). Then, alignments were trimmed for the conserved domain manually on Geneious and automatedly trimmed using trimAl (95, 96). The tree was made using RAxML (97) with the JTT protein substitution model (98). The trees were analyzed and figures were generated using iTOL (99). The 18s genes were found and the alignment was generated using ssu_tree.py (https://github.com/christophertbrown/bioscripts27), which searches for rRNA genes using an HMM method then the sequences were dereplicated and clustered at 98% nucleic acid identity, again representing a possible species level designation (100, 101), and aligned using SSU-ALIGN (102). Assembly errors were found and fixed in the scaffolds containing the 18S using ra2.py (https://github.com/christophertbrown/fix_assembly_errors) (90). The tree was generated using the approach described above.

We identified the RNA viral scaffolds using a combined profile HMM approach, using HMMs from Pfam (86) for different types of RdRps. The Leviviridae were identified using RNA_replicase_B (PF03431) and scaffolding errors were identified with ra2.py. There were no useable HMMs for the Cystoviridae, so we generated our own. We aligned public RdRp sequences from the Cystoviridae using MAFFT v7.402 (E-INS-i) on Cipres and generated an HMM using hmmbuild from the HMMER suite with default settings (87, 93, 94). The Cystoviridae HMM is publicly available on the figshare repository corresponding to this report. We used many publicly available HMMs to find the RdRp of eukaryotic RNA viruses (Mononeg_RNA_pol [PF00946], RdRP_5 [PF07925], Flavi_NS5 [PF00972], Bunya_RdRp [PF04196], Mitovir_RNA_pol [PF05919], RdRP_1 [PF00680], RdRP_2 [PF00978], RdRP_3 [PF00998], RdRP_4 [PF02123], RVT_1 [PF00078], RVT_2 [PF07727], Viral_RdRp_C [PF17501], and Birna_RdRp [PF04197]). The RdRp genes were initially classified using a blast and tree-building method to determine the correct translation table to use for Prodigal in the single mode (85). For mitoviruses we used translation table 4, which to our knowledge is used by all mitoviruses. However, the additional mitoviral genes we predicted may have been incorrectly called if these genomes use a modified genetic code, as occurs in some fungi (103, 104). We examined gene predictions for indications of this (e.g., interruption of the RdRp gene) but as this phenomenon was not identified, no alternative codes were used for gene prediction. RdRp amino acid sequences were dereplicated and clustered at 99% amino acid identity using USEARCH (92). We used previously published alignments for many viral families (36) and added our sequences and key references using the Mafft v7.402 with the seed (93) and E-INS-i options on Cipres (94). For viral clades without published alignments or where the published alignment was inappropriate (Fusariviridae, Narnaviridae, Leviviridae, and Cystoviridae) we generated our own alignments using reference sequences from NCBI and the same alignment and tree-building steps as described above.

To identify lysis proteins in the Leviviridae genomes we used the Geneious ORF prediction to find all possible ORFs (96). The amino acid sequences were run through PSORTb v3.0.2 with the gram-negative setting to find possible lysis proteins (105).

To obtain coverage values for the viruses, we mapped against the entire scaffold containing the RdRp gene. For the presumed hosts we only mapped reads to the ORF, for the rpS3 and Cox1. Raw reads from all samples were mapped using Bowtie2 (–sensitive and–rfg 200,300 options), then reads were filtered for 2 mismatches using calculate_breadth.py (https://github.com/banfieldlab/mattolm-public-scripts/blob/master/calculate_breadth.py) (106). Coverage values were converted to read counts and normalized with DESeq2 (107). Ordinations were generated from DESeq2 normalized count data in R. The data were ordinated using nonmetric multidimensional scaling (R package: vegan) and significantly different clusters were determined using adonis (108). Using DESeq2 we determined the viral and eukaryotic habitat enrichments in the treatments (rhizosphere+litter, bulk+litter, and rhizosphere) relative to bulk soil at each time point. We also conducted the opposite test, identifying enrichments in bulk soil by comparing normalized abundance to each treatment (rhizosphere+litter, bulk+litter, and rhizosphere) at each time point. P values were corrected for multiple comparisons using a Benjamini–Hochberg correction (109).

Network analysis was conducted using the DESeq2 normalized coverage for the eukaryotic viruses and bacterial viruses and the eukaryotic Cox1 and bacterial rpS3. These values were filtered for coverage greater than 1,000. Spearman’s rank correlation coefficient values were computed in R using the Hmisc package; the P values were corrected for multiple testing (109, 110). The data were then stringently filtered to contain only organisms present in 20 or greater samples, a corrected P value of 0.0001 or less, a correlation value of 0.7 or greater, and removed all host–host correlations. These data were then fed into Cytoscape and a network was constructed using default settings (111). Only clusters that contained at least 1 virus as the dependent variable and 1 host as the independent variable were analyzed and presented.