ABSTRACT
Biological soil crusts (biocrusts) are photosynthetic “hot spots” in deserts and cover ∼12% of the Earth’s terrestrial surface, and yet they face an uncertain future given expected shifts in rainfall events. Laboratory wetting of biocrust communities is known to cause a bloom of Firmicutes which rapidly become dominant community members within 2 days after emerging from a sporulated state. We hypothesized that their bacteriophages (phages) would respond to such a dramatic increase in their host’s abundance. In our experiment, wetting caused Firmicutes to bloom and triggered a significant depletion of cyanobacterial diversity. We used genome-resolved metagenomics to link phage to their hosts and found that the bloom of the genus Bacillus correlated with a dramatic increase in the number of Caudovirales phages targeting these diverse spore-formers (r = 0.762). After 2 days, we observed dramatic reductions in the relative abundances of Bacillus, while the number of Bacillus phages continued to increase, suggestive of a predator-prey relationship. We found predicted auxiliary metabolic genes (AMGs) associated with sporulation in several Caudovirales genomes, suggesting that phages may influence and even benefit from sporulation dynamics in biocrusts. Prophage elements and CRISPR-Cas repeats in Firmicutes metagenome-assembled genomes (MAGs) provide evidence of recent infection events by phages, which were corroborated by mapping viral contigs to their host MAGs. Combined, these findings suggest that the blooming Firmicutes become primary targets for biocrust Caudovirales phages, consistent with the classical “kill-the-winner” hypothesis.
IMPORTANCE This work forms part of an overarching research theme studying the effects of a changing climate on biological soil crust (biocrust) in the Southwestern United States. To our knowledge, this study was the first to characterize bacteriophages in biocrust and offers a view into the ecology of phages in response to a laboratory wetting experiment. The phages identified here represent lineages of Caudovirales, and we found that the dynamics of their interactions with their Firmicutes hosts explain the collapse of a bacterial bloom that was induced by wetting. Moreover, we show that phages carried host-altering metabolic genes and found evidence of proviral infection and CRISPR-Cas repeats within host genomes. Our results suggest that phages exert controls on population density by lysing dominant bacterial hosts and that they further impact biocrust by acquiring host genes for sporulation. Future research should explore how dominant these phages are in other biocrust communities and quantify how much the control and lysis of blooming populations contributes to nutrient cycling in biocrusts.
INTRODUCTION
The Southwestern United States is predicted to experience alterations in the timing of rainfall events with concomitant increases in soil surface temperatures (1). In arid lands, moisture pulses are infrequent, and yet their timing, intensity, and duration dictate biological activity (2). A significant avenue of recent microbial ecological research has been elucidation of the responses of microbial communities to various manifestations of climate change, such as increased rainfall, and of how these responses in turn affect ecosystem processes (3–5).
Biological soil crusts (referred to as biocrusts) are dominant topsoil microbial communities in deserts and serve as the basis for soil surface colonization, stabilization, and activity (6). Biocrusts can be comprised of photosynthetic cyanobacteria, bryophytes (mosses and liverworts), lichens, fungi, a myriad of heterotrophic lineages that cooperate to cycle both soil carbon and nitrogen, and other elements (7). Our previous research in biocrusts dominated by cyanobacteria has shown that biocrust wetting results in significant shifts in microbial community structure (8) and dramatically alters microbial metabolic outputs within the first day of wetting (9). The most dramatic of these observations is that multiple Firmicutes lineages bloom and then collapse following wetting, which we hypothesize is due to bacteriophage (phage) predation.
Recent studies using cultivation-independent metagenomic sequencing have revealed the roles that viruses play across terrestrial biomes, including predation by giant viruses in forests (10) and viral predation of dominant and ecologically relevant microbial lineages along a permafrost thaw gradient (11, 12) and of dominant tailed phages in Antarctic soils (13, 14). Viral metagenomics offers insight into viral diversity within complex samples (15) and unlocks the potential for novel viral diversity to be uncovered without the need for cultivation (16).
We have previously shown that Firmicutes spp. respond positively to a biocrust wetting event (8), which is thought to be due to their emergence from dormancy (17). However, we noted previously that many Bacillus strains declined in abundance soon after the initial wetting event whereas others continued to increase in abundance until a sudden decline occurred (9). We hypothesized that the increase in phage activity upon wetting leads to the predation of blooming Firmicutes populations in biocrust. Several lines of evidence motivate this view: First, phages that infect Bacillus in hot desert soils in both the Mojave Desert and the Namib Desert have been described previously (18, 19). Second, phages have been shown to target blooming bacterial populations in aquatic ecosystems (20), including Prochlorococcus spp. in the ocean (21, 22) and blooms of the cyanobacterium Microcystis aeruginosa in freshwater lakes (23). Finally, the “kill-the-winner” hypothesis posits that viruses act as moderating agents within a community by attacking the most active population of prokaryotes, thus making the most of the limited resources available to them (24, 25). Accordingly, we expected biocrust phages to target the most active hosts following a latent period. To test this concept, we mined microbial and phage sequences from metagenomes obtained from biocrust exposed to a controlled laboratory wetting event and tracked their responses over 2 days. The samples originated from the Moab Desert, a site of long-term research aimed at exploring the effects of climate change in arid regions in Utah in the Southwestern United States.
RESULTS
Our 24 biocrust communities were subjected to a simulated wetting experiment, which we leveraged to explore microbial and phage responses along the 2-day time series. The assembled sequence data comprised 28.4 million contigs (∼20 Gbp) (see Table S1A and B in the supplemental material) from which we reconstructed 40 nearly complete bacterial metagenome-assembled genomes (MAGs). According to the Genomic Standards Consortium standards (26), all of the MAGs were of medium quality, as most lacked an identifiable 16S rRNA gene, and were >75% complete with <10% contamination (Table S2). These genomes encompassed bacterial members of Firmicutes (n = 16), Bacteroidetes (n = 7), Proteobacteria (n = 5), Acidobacteria (n = 4), Chloroflexi (n = 2), and Actinobacteria (n = 2) and included a single representative of Deinococcus-Thermus (see Fig. S1A and B in the supplemental material). The assembled genomes averaged 4.4 Mbp in size, with a range of relative abundances from high abundance (1.33% of assembled metagenome contigs) to rare (0.0067%) community members. We interrogated our MAGs for genes of interest relating to nutrient cycling, sporulation, osmoregulation, carbon starvation, and heat shock tolerances (see Results in Text S1 in the supplemental material).
TEXT S1
This is a work of the U.S. Government and is not subject to copyright protection in the United States. Foreign copyrights may apply.
FIG S1
This is a work of the U.S. Government and is not subject to copyright protection in the United States. Foreign copyrights may apply.
TABLE S1
This is a work of the U.S. Government and is not subject to copyright protection in the United States. Foreign copyrights may apply.
TABLE S2
This is a work of the U.S. Government and is not subject to copyright protection in the United States. Foreign copyrights may apply.
The most phylogenetically diverse biocrust communities were the prewetting (“dry”) metagenomes and those collected 3 min and 9 h after wetting (∼58% of the community was sequenced on average as estimated using Nonpareil [27]). Biocrust metagenomes from later time points, 18 and 42 h after wetting, captured ∼64% of the community, while the final time point, 49.5 h after wetting, produced communities that were less complex, whereby 71% of the community was sequenced (Fig. S2A). Diversity estimates corroborated these rarefaction estimates and showed significant decreases in α-diversity following wetting (analysis of variance [ANOVA], P < 0.05) (Table 1). The 49.5-h group had significantly lower α-diversity than any of the other groups (Tukey post hoc test, P < 0.05). In contrast, the number of single-copy ribosomal genes (large and small subunits) increased significantly with wetting (P < 0.05). We also observed a significant increase in β-diversity over the course of our wetting experiment (within treatment variation). Note that the decrease in diversity may also have been due to our limited sampling and sequencing of rare community members as their relative abundances decreased dramatically after wetting.
Statistics describing the taxonomic diversity within the biocrust metagenomes grouped according to time after wetting
FIG S2
This is a work of the U.S. Government and is not subject to copyright protection in the United States. Foreign copyrights may apply.
We observed significant differences in community composition according to time after wetting (permutational multivariate analysis of variance [PERMANOVA], P < 0.05) (Fig. 1, left panel). The effect of wetting induces changes in composition that are strong enough to supersede those of the maturity stage of the biocrust, which developed over decades (P > 0.05). The increase in β-diversity over time was also observable as individual communities separated within ordination space (nonmetric multidimensional scaling [NMDS] ordinations based on Bray-Curtis dissimilarities of 16S rRNA gene content) (Fig. 1, left panel), which was the result of taxonomic divergence over time and of overall increases in 16S rRNA gene content (Fig. 1, right panel).
(Left) The changes in taxonomy over the course of the wetting experiment are observable in NMDS ordination space. Biocrust communities grouped according to time after wetting, which was sufficient to explain the differences in community composition (PERMANOVA; P < 0.05). (Right) The communities grew more dissimilar as the copy numbers of 16S rRNA genes increased over time.
Overall, the biocrust communities encompassed 23 bacterial phyla and 1 archaeal phylum (Thaumarchaeota). The metagenomes also provided contigs belonging to Ascomycota (filamentous, saprobic Fusarium and Purpureocillium lineages), metazoa (Tardigrada), and mosses (Physcomitrella patens). The low diversity of fungi might have been due to database bias, and using a dedicated fungal database such as UNITE might have revealed more fungal diversity in these samples. Next, we explored the taxonomic shifts that occurred during biocrust wetting by tracking changes in 16S rRNA gene content in the unassembled reads. On the basis of our previous 16S rRNA gene survey, we anticipated a significant enrichment of Bacillales and other Firmicutes taxa ∼18 h after wetting (8, 9). Notably, these blooms should have been accompanied by substantial reductions in the proportion of filamentous cyanobacterial community members across all samples as nonbloomers are outcompeted by faster growers. However, Cyanobacteria are known to be sensitive to wetting (28, 29), and excessive cellular swelling and lysis are thought to cause their negative response. The dry and 3-min biocrust metagenomes were dominated by photosynthetic Cyanobacteria and by heterotrophic Actinobacteria, Bacteroidetes, and Proteobacteria (Fig. 2, left panel). Overall, the biocrust communities encompassed 23 bacterial phyla and 1 archaeal phylum (Thaumarchaeota). The metagenomes also provided contigs belonging to Ascomycota (filamentous, saprobic Fusarium and Purpureocillium lineages), metazoa (Tardigrada), and mosses (Physcomitrella patens).
Exposure to a sudden wetting event had strong effects on the taxonomy of the biocrust communities across all successional stages. (Left) Relative abundances of microbial phyla in the biocrust communities. Archaea and eukaryotes are classified at the domain level, while Proteobacteria are distinguished at the class level. (Right) The responses of Firmicutes genera are shown as the proportion of 16S rRNA gene copies per metagenome.
The wetting event drove decreases in the relative abundances of all nonblooming phyla. Notably, the relative abundance of Cyanobacteria decreased significantly from ∼64% of the community in the dry samples to just 0.23% 49.5 h after wetting (P < 2.30e−10). We observed striking, significant decreases in the relative abundances of most cyanobacterial lineages over the course of the experiment, culminating in a nearly complete loss of the genera Microcoleus (P < 2.15e−7), unclassified Cyanobacteria (P < 4.78e−5), Tychonema (P < 2.36e−3), and Subsection III, Family I members (P < 0.025) by 49.5 h (Fig. 2, left panel). Cyanobacteria are known to be sensitive to wetting (28, 29) and excessive cellular swelling, and lysis may have caused their negative response here. The reduction in the abundances of cyanobacterial members also had dramatic effects on network topology across the experiment on the basis of co-occurrence analyses of taxonomic information (see Text S1, Fig. S3A to C, and Table S5).
FIG S3
This is a work of the U.S. Government and is not subject to copyright protection in the United States. Foreign copyrights may apply.
TABLE S5
This is a work of the U.S. Government and is not subject to copyright protection in the United States. Foreign copyrights may apply.
Only members of the Firmicutes significantly increased in relative abundance after wetting; their levels increased from 1.07% of 16S rRNA genes in dry samples to 82.3% of 16S rRNA genes at 49.5 h (relative abundance of 16S rRNA genes, ANOVA, Bonferroni corrected, P < 1.77e−6) (Table S3A). The presence of the Firmicutes bloom was pronounced across all biocrust successional stages, but the responses of individual genera were more nuanced. Notably, Brevibacillus showed a steady increase in abundance after wetting (Fig. 2, right panel), whereas Paenibacillus showed slower growth. The population of Bacillus became dominant 9 h after wetting but had collapsed by 49.5 h. We speculated that these diverse responses to wetting would be reflected by increases in the proportions of functional genes affiliated with the changes in the proportions of these taxa.
TABLE S3
This is a work of the U.S. Government and is not subject to copyright protection in the United States. Foreign copyrights may apply.
We tracked how changes in taxonomic composition affected the relative contributions of functional genes. Most of the upregulated genes in our analysis were involved in sporulation (Fig. 3), suggesting that previously dormant Firmicutes had replicated within 9 h after wetting. Specifically, we observed significant increases in the abundance of sporulation cluster stage IIIA genes (AB, AC, AD, and D), stage V genes (spoVS), spore maturation protein A, and genes encoding the spore cortex-lytic enzyme CwlJ (Fig. 3). Consistently, the dramatic reduction in the populations of cyanobacterial taxa over the course of the wetting experiment led to significant reductions in the relative abundances of genes involved in the assembly of photosystem I (ycf4 and psaLK2) and photosystem II (psbADMPVY) and of genes involved in the production of bilin (pcyA) and the phycobilisome (nblB and cpeS). The relative abundances of genes associated with the cyanobacterial circadian clock (input kinase A) and myxoxanthophyll biosynthesis (cruF) and genes encoding products associated with photorespiration (d-glycerate 3-kinase and glycolate oxidase) were also significantly reduced (Table S3B). Strikingly, numerous genes involved in carbon dioxide (CO2) uptake (carboxysome) were also significantly reduced in abundance. These included RuBisCO chaperonin genes (rbcX) as well as CO2-concentrating genes (ccmO), high-affinity carbon uptake genes (hatR), and low-affinity CO2 hydration genes (cphX). Overall, 95 genes decreased in relative abundance over the wetting experiment, while only 22 genes became significantly more abundant following wetting (Table S3B).
Heat map indicating the distribution of major functional genes in biocrust metagenomes following soil wetting. Count data are log normalized to enable comparisons across different sequencing efforts. Genes that are overrepresented are represented by darker green colors; white data represent genes that are underrepresented (i.e., the genes with the lowest relative abundance levels). The x axis indicates biocrust communities following wetting, i.e., ranging from the dry state to 49.5 h after wetting. The y axis data include all functional genes that differed significantly in abundance across the wetting event.
Next, we used VirSorter to mine viral signals from our 28 metagenomes (including the 4 bundle metagenomes), and the mined signals were coassembled into 5 data sets to assist viral sequence recovery. This procedure produced a total of 1,849 unambiguous viral sequences (see Materials and Methods). The removal of short (<10-kb) contigs, together with manual curation (see Materials and Methods), refined our data set to 1,210 viral genomes (ViGs) or genome fragments.
To place our ViGs within the known, relevant diversity of publicly available viral genomes, we created a virus genome-based network based on shared protein content (30, 31) by combining our contigs with the viral genomes and genomic fragments present in the RefSeq viral database (v85) (32). The analysis produced 40,796 virus clusters (VCs) representing approximately genus-level taxonomy, of which 2,294 (6%) contained biocrust viral contigs, representing 504 uncultivated viral genomes (UViGs). Overall, 148 (29%) of the 504 biocrust UViGs that we identified were found to be taxonomically affiliated using these VCs. Clustering the 504 UViGs at 95% average nucleotide identity (ANI) and 85% alignment fraction (AF) (20) produced 433 viral operational taxonomic units (vOTUs), of which 49 could be assigned at the genus level. All classified vOTUs belonged to Caudovirales families (double-stranded DNA [dsDNA] tailed phages), although six vOTUs were broadly categorized as representing dsDNA viruses (Fig. 4, left and top right panels; see also Table S4A).
(Left) Network clustering of recovered virus genomes and genome fragments with RefSeq viral genomes (v85) based on shared predicted protein content. Major viral families (nodes) are indicated by unique colors and shapes, with vOTUs found in this study indicated by green circles. Edges (lines between nodes) indicate significant shared protein content. (Top right) Stacked bar chart indicating the number of UViGs that were recovered across the wetting experiment. (Bottom right) Correlation between the relative proportions of Firmicutes and UViGs per time point.
TABLE S4
This is a work of the U.S. Government and is not subject to copyright protection in the United States. Foreign copyrights may apply.
The longest circular UViGs belonged to three identical Spbetavirus phages (163,801 bp), which are Siphoviridae phages known to infect soil Bacillus (33). These UViGs were reconstructed independently biocrust metagenomes seen in the early, late-middle (late-mid), and late successional stages, while the early-middle (early-mid) metagenome also contained a smaller Spbetavirus contig of 137,161 bp. the Siphoviridae phages were the most abundant, with 89 vOTUs spanning 17 genera. Podoviridae were represented by 19 vOTUs, including three T7 viruses, while 21 vOTUs assigned to the Myoviridae family were found to include three genera (Abouovirus, Bxz1virus, and P2virus). A further 10 genomes were coarsely assigned to the Caudovirales. Mapping UViGs back to individual metagenomes revealed that all of the phages were detected only after wetting; i.e., no phages were recovered from the bundles or the dry samples (Fig. 4, top right panel), which culminated in the greatest number of UViGs (n = 50) being recovered more than 2 days after wetting. For example, the full-length spbetaviruses emerged at 42 h and 49.5 h after wetting, which correlates with the decline in the abundances of the putative Bacillaceae hosts at those time points (Fig. 2, right panel). Andromedaviruses (length range, 47,888 to 48,100 bp) were found in all successional stages at 18, 42, and 49.5 h after wetting and have previously been shown to infect Bacillus hosts (33). Combined, these UViG data indicated that Firmicutes were predominant targets, particularly at the end of the wetting experiment, which led us to further explore this putative link. The classification of our viral genomes provided primary evidence of phages targeting Firmicutes such as Bacillus, Staphylococcus, and Brevibacillus.
Putative auxiliary metabolic genes (AMGs; genes that originate from microbes and can increase fitness) that code for sporulation proteins (spoIIIE) were identified in 9 UViGs. To confirm that the AMGs were of viral origin, we checked that flanking regions present both upstream and downstream contained common virus genes (integrases, structural genes, or terminases), had significant matches to the gene of interest (E value < 1 × 10−3; bit score > 60), and represented the only AMG per contig. BLASTP queries indicated that all 9 sporulation genes (spoIIIE) were assigned to Bacillus strains. We used the sporulation AMGs to explore the proportion of phages infecting the blooming Firmicutes. As an initial proof of concept, we found that the numbers of UViGs correlated strongly with the relative abundances of Firmicutes across the wetting experiment (r = 0.762). That is, as the proportion of Firmicutes increased throughout the wet-up, so too did the number of recoverable phage genomes (Fig. 4, bottom right panel; see also Fig. S4B).
FIG S4
This is a work of the U.S. Government and is not subject to copyright protection in the United States. Foreign copyrights may apply.
Genes for sporulation were detected in 15 of the 16 Firmicutes MAGs. Genes involved in spore DNA protection (dipicolinate synthase subunits A and B, dpaAB) and small acid-soluble spore proteins (SASP1, SASPS6, and gpr) were detected exclusively in the Firmicutes genomes and are well-known components of endospore production (34, 35). DpaAB genes code for dipicolinic acid, which allows endospores to maintain dormancy under harsh conditions. Sporulation-related genes, including spoIIIE and gerE, were among the most highly overrepresented genes in the metagenomes examined in our study by the end of our wetting experiment (Tables S3B and S4B), reiterating the increase in the relative abundances of spore-formers following wetting. A maximum likelihood tree indicated that the spoIIIE AMGs found in our UViGs showed high similarity to those found in the Firmicutes genomes (Fig. 5). We confirmed the clustering result by predicting protein structures of each spoIIIE gene using I-TASSER and comparing models using TM-align, which showed a similarity score of 0.488.
Maximum likelihood tree indicating the phylogenetic placement of spoIIIE genes retrieved from UViGs with their nearest bacterial neighbors (from Firmicutes MAGs). Protein structures were predicted using I-TASSER, and the first predicted models are presented. Tree lineages are colored blue for bacterial sequences and orange for viral sequences.
Additionally, we searched for CRISPR-Cas repeats and proviral genes within the reconstructed MAGs to link viruses to potential hosts. We found evidence of CRISPR-Cas repeats in 16 MAGs, which links them to previous or ongoing infection events by biocrust viruses (11, 36). Six of the 16 Firmicutes MAGs contained CRISPR elements (Table S2). Aligning the CRISPR-Cas repeats and spacers from the MAGs with our virus contigs revealed perfect homology between a Phietavirus (EM_51679) and the Anoxybacillus MAG EM_40. Both the repeat sequence (5′-CTTCACATCACACATAGTTTCCTCAGAAAC-3′) and three spacers (data not shown) were identical between the virus and its cognate host, which offers a direct link to the microbial host in silico. Next, we searched for prophage elements within each MAG to uncover temperate phages that utilized a lysogenic infection cycle. We found 12 MAGs with viral elements within their genomes, of which 7 were Firmicutes (Table S2).
DISCUSSION
The results presented here provide an initial assessment of the phage diversity in biocrust communities—dominant photosynthetic topsoil communities found around the globe. This description of dominant phage lineages in biocrust further extends our knowledge of viruses in important soil ecosystems (37), which lags considerably behind what is known about marine viruses (20, 38). Our work sought to resolve the contradicting fates of taxa within the Firmicutes after biocrust communities were exposed to a simulated wetting event (8), which we predicted might have been due to phage predation. Here, we provide evidence of a predator-prey relationship between Bacillus and their phages. Comparing the 16S rRNA gene copy numbers of Bacillus throughout wet-up to the numbers of their known UViGs showed that phage abundance peaked as that of their hosts began to decline (Fig. 6).
Predator-prey relationship between Bacillus organisms and their phages. Bacillus rRNA copy number is provided on the left-hand y axis with error bars calculated across four metagenomes from each time point. The numbers of UViGs linked to Bacillus infection are indicated on the right-hand y axis.
The presence of proviral sequences within Firmicutes genomes suggests that lysogenic replication strategies are common among biocrust viruses. Temperate viruses shift their replication strategy from lysogenic to lytic replication as bacterial host production increases (39). We found prophages in 12 of the 40 MAGs. Seven Paenibacillus MAGs contained proviral sequences, and each of these contained multiple prophages. While the genus Bacillus abundance declined, both Paenibacillus and Brevibacillus continued to increase in relative abundance throughout the wet-up. According to the kill-the-winner model, we would expect phages to attack the fastest-growing host population first (i.e., Bacillus). We anticipate that phages would target members of the Paenibacillus and Brevibacillus after the Bacillus populations have been depleted on the basis of the trajectory of these genera over 2 days (Fig. 2, right panel).
For Bacillus hosts, phage-encoded sporulation genes could play an additional role in these dynamics by providing a mechanism for lysis with concomitant phage and spore release where uninfected host spores would represent a reservoir of future host cells. Note that we cannot differentiate between population reduction mediated by lysis from sporulation and that mediated by some combination of factors with these data. In addition, spores that include integrated phage genomes would also provide additional protection, which may be important in the harsh biocrust environment. Regardless, the lysis of Firmicutes cells probably stimulates the biocrust community as metabolites, such as nitrate, are released into the soil. Conceivably, this viral shunt could provide essential nutrients to mosses or other higher eukaryotes, as predicted in marine ecosystems (40), which has the potential to influence ecological services of the crust.
Our metagenome-assembled genomes also revealed that Firmicutes-infecting phages acquired genes horizontally from their hosts’ genomes. Specifically, nine spoIIIE genes in phage genomes link them to Firmicutes infections. The spoIIIE gene codes for a DNA translocase that plays a critical role during sporulation by assisting the translocation of chromosomal DNA from mother cells into the forespore during the final stages of engulfment (41). It has been shown that marine cyanophages contain photosynthetic genes which maintain essential host functions (22). These spoIIIE genes may have a similar role in these biocrusts if sporulation is an important component of the life history of these Bacillus phages, as endospores can be carriers of phage genomes, thus offering protection from various environmental stresses (42, 43). Alternatively, it is possible that these phage-encoded proteins were originally sporulation related but now perform another function related to DNA export and/or ATP-dependent motor activity. AMGs encoding sporulation machinery (spoVS) in vOTUs from permafrost have been detected previously (12), and it would be intriguing to investigate how widespread these genes are in other soil systems as they may influence host sporulation dynamics (42).
As Firmicutes bloom during precipitation events, they are probably preyed upon by phages, exerting control on their abundance throughout the year. Biocrust saturation after rainfall releases adsorbed phages, which may initiate repeated cycles of lytic replication (37). The drying soil likely draws phages into small pores by shrinking water films (44). Soil microheterogeneity may further limit lytic virus-host interactions as cells become heterogeneously distributed within the dry soil matrix. This paradigm is especially true of desert biocrusts, which experience prolonged periods of dormancy for much of the year and reawaken after moisture events (45, 46). Future research should investigate these dynamics in the field, particularly as our long-term storage of samples (∼11 months in the dark) may have influenced the stress levels that the microbes experienced. Moreover, it is unclear whether blooms of Firmicutes occur as rapidly in naturally occurring biocrust communities or if this response is mitigated by other factors.
Previous surveys of viral populations from desert biotopes, such as cyanobacterium-dominated hypoliths (18) and exposed desert soils (47), revealed a high proportion of Caudovirales (dsDNA tailed phages), which appear to be dominant dsDNA phages in biocrust as well. We found that the family-level hierarchical abundance of Siphoviridae > Myoviridae > Podoviridae was remarkably similar to that described in Antarctic hypoliths (14), potentially indicating that these tailed phages are core constituents of “hot spots” of microbial activity in desert ecosystems. Moreover, the average genome sizes of phages identified here (average length, ∼43 kb) are like those described in Namib soils identified from pulsed-field gel electrophoresis profiles (55 kb). Additionally, those studies provided striking evidence that cyanophages are a potentially low-abundance component of desert viral communities despite the abundance of putative cyanobacterial hosts in hypoliths (47) and biocrust, including the genera Microcoleus, Scytonema, Schizothrix and Nostoc (7, 48). These results should be interpreted with caution, however, as our long-term storage of biocrust samples in the dark for 11 months may have diminished the responsiveness of certain taxa to wetting. Moreover, most knowledge regarding cyanophage diversity stems from oceanic surveys, which may have introduced biases into the RefSeq database that was used here. The number of studies on soil viruses continues to increase rapidly, which certainly offers researchers the potential to detect novel cyanophages in these underexplored soil ecosystems.
Summary.Biocrust wetting induced a bloom of Firmicutes and their phages. Our data suggest that Firmicutes emerge from spores, whereby fewer individuals remain in sporulated states after 2 days, which is reflected in significant enrichment for genes relating to sporulation. In contrast, the reduction in the relative abundance of Cyanobacteria resulted in significant depletion of functions relating to both carbon and nitrogen fixation, and this was observed in collapses in co-occurrence network structures (see Fig. S3C in the supplemental material). Our investigation of phage diversity in biocrust, performed here for the first time to our knowledge, did not link phage predation to cyanobacterial decline but instead revealed that biocrust phages primarily infect blooming Firmicutes, particularly members of the Bacillus genus. The phage-bacterium dynamics observed here appear to follow the kill-the-winner model, in which phages expand on the fastest-growing host population in biocrust. The collapses of population blooms are often mediated by phage infection and lysis and can occur on very short time scales such as in marine systems where phages target blooming Synechococcus species within 30 min (40). Here, we showed that the increase in Firmicutes abundance over 2 days is closely followed by peaks in the frequency of recoverable Caudovirales phage genomes targeting different members of the phylum, followed by a decline of specific Firmicutes lineages as their phages replicate. The presence of putative auxiliary metabolic genes associated with sporulation in these phages suggests that the sporulation process may be an important component of the life history of these phages. Future research could further investigate the contributions of sporulation to phage dispersal and maintenance within endospores. Additional analyses could consolidate the fate of the biocrust communities by investigating how biocrusts return to their initial steady state following drying.
MATERIALS AND METHODS
Sample collection, experimental setup, and DNA extraction.Details of the sample collection procedure, the biocrust wetting experiment, and DNA extraction protocols have been described previously (9) and are illustrated in Fig. S4A in the supplemental material. Briefly, 28 petri dishes (6 cm2 by 1 cm in depth) were used to core biocrust samples from the Green Butte Site near Canyonlands National Park (38°42′54.1′′N, 109°41′27.0′W, Moab, UT, USA) in 2014. The site is an area of long-term ecological research interest aimed at exploring climatic changes in arid regions. Biocrust samples (28 petri dishes) were collected along a maturity gradient of Cyanobacteria-dominated biocrusts ranging from lightly colored, young crust to dark mature crust, which encompassed four maturity stages overall. Dry samples collected in the field were maintained in a dark desiccation chamber in the laboratory for 11 months until required. Biocrust samples underwent a laboratory (12 h light/12 h dark cycle) wetting event by wetting of each soil sample (0.5 g) with liquid chromatography-mass spectrometry (LC-MS)-grade water (1 ml). Soil was removed at five time points following wetting (in addition to dry and bundled samples) for each successional stage (3 min, 9 h, 18 h, 42 h, and 49.5 h after wetting). After pelleting of soil by centrifugation (5,000 × g for 5 min) and removal of porewater, DNA was extracted from soils (0.25 g) using a DNeasy PowerSoil kit (Qiagen, Hilden, Germany), resulting in 100 μl of eluted DNA. Library preparation and sequencing were performed at the QB3 facility at the University of California, Berkeley (UC Berkeley), using an Illumina HiSeq 4000 system.
Metagenomic sequencing and analysis.Metagenomic sequencing provided ∼413 Gbp of Illumina HiSeq data in 2,754,282,840 raw sequences across 28 biocrust samples. The data included four separate metagenomes for microbial communities collected before wetting and those collected 3 min, 9 h, 18 h, 42 h, and 49.5 h after wetting. Sequencing was performed on an Illumina HiSeq 4000 instrument using 2 × 150-bp chemistry. Raw sequences were filtered using Prinseq v0.20.4 (49), which eliminated reads with phred scores of <20 for six consecutive bases (–min_qual_mean), singletons, and those with ambiguous bases (ns_max_n). Estimates of sequence coverage were obtained using Nonpareil v3.30 (27). High-quality sequences were assembled under default conditions using MEGAHIT v1.1.3 (50) as previously suggested for highly complex soil samples (51). These 24 metagenomes were also coassembled according to their maturity stage, i.e., dry, early, early-mid, late-mid, and late successional stage biocrust, as physically identified a priori during the field sampling expedition. Statistics about each assembled metagenome were obtained using MetaQUAST v4.6.3 (52) and are available in Table S1B in the supplemental material. Open reading frames (ORFs) were extracted from each assembled metagenomes with Prodigal v2.6.3 (53). Annotation of contigs was performed using DIAMOND v0.9.14.155 (54) with queries made against the entire NCBI nonredundant nucleotide database. Here, annotations were provided at E values of <0.00001 (-e), while the five best hits were retained (-k). MEGAN software v6.12.5 was used to visualize functional and taxonomic attributes of communities (55). We extracted all 16S rRNA genes from the curated metagenomes using SortMeRNA (56). Finally, we resolved exact sequence variants (ESVs) using DADA2 (Divisive Amplicon Denoising Algorithm [57]), with default parameters except for truncLen (120, 150) and maxEE (2, 5). We assigned taxonomy against the entire SILVA reference database (58) after random rarefaction of each community to the lowest sequencing effort (8,805).
Network construction.We constructed co-occurrence networks at the genus level using our log-transformed taxonomic data. Networks were constructed using Sparse Correlations for Compositional (SparCC) data (59). We chose to divide our 24 metagenomes according to time after wetting. Individual networks were created for the beginning (dry and 3 min after wetting), middle (9 h and 18 h after wetting), and end (42 h and 49.5 h after wetting) of the experiment. We used significant SparCC values (P < 0.05) that were greater than 0.6 or less than −0.6 and included only genera present in at least half of the samples within our OTU table. All networks were visualized using Cytoscape v3.6.1 (60), and highly connected modules were automatically identified using MCODE (61).
Reconstructing genomes from metagenomes.Due to our ultradeep sequencing effort, we were able to reconstruct metagenome-assembled genomes (MAGs) from our assembled biocrust communities using MetaBAT v2.12.1 (62). We increased the sensitivity of the binning process by increasing the minimum contig length to 5,000 bp (–m), the percentage of good contigs to 97 (–maxP), the minimum edge binning score to 85 (–minS), and the maximum number of edges per node to 250 (–maxEdges). We also included the –noAdd option, which prevented small contigs from being binned and prevented the creation of bins that were <1 Mbp in size. MAGs were checked for completeness and contamination using CheckM v1.0.8 (63). Medium- to high-quality MAGs were retained according to the conventional values for minimum information about a MAG (>75% completeness, <10% contamination) (26). MAGs were then functionally annotated online using Prokka v1.12 (64) and the RAST v2.0 server (65), with ORFs validated with UGENE (66). Quality statistics about the MAGs are presented in Table S2.
Phylogenetic tree building.Phylogenetic trees were built using universal bacterial protein gyrase subunit B (gyrB) extracted from MAGs (n = 19) and 2,878 reference sequences. The gyrB sequences were collated and aligned using MAFFT v7 (67). Neighbor-joining trees were built with bootstrap values reported, and visualization was performed in iTOL v3 (68). The phylogenetic placement of the MAGs was also evaluated using the Microbial Genome Atlas (MiGA) online (69).
Statistical analyses.We used STAMP v2.1.3 (70) and R v3.5.0 to make statistical comparisons between the communities according to both time after wetting and successional stage. Within the R statistical environment, we used the vegan package (71) to explore ecologically meaningful differences between communities according to time after wetting and successional stage. Bray-Curtis dissimilarity matrices were used to explore biocrust community structure after Hellinger transformations of species relative abundance data. PERMANOVA was used to test for differences in community structure on the basis of time after wetting (dry biocrust and biocrust at 3 min, 9 h, 18 h, 42 h, and 49.5 h after wetting) and successional maturity (early, early-mid, late-mid, and late) using the adonis function. We used STAMP to test for significant differences in taxonomic and functional composition between time points and successional stages and used ANOVA with Tukey-Kramer post hoc tests and Bonferroni corrections of the P value statistic.
To obtain viral contigs from biocrust metagenomes, we coassembled samples according to their characteristics and community complexity. The coassemblies were produced by combining the 4 filamentous bundles, the 4 dry biocrust samples, the 5 early successional stage biocrust samples, the 5 early-mid succession biocrust samples, the 5 late-mid succession biocrust samples, and the 5 late succession biocrust metagenomes into single, concatenated metagenomes. This produced 6 coassemblies to complement the 28 individual metagenomes.
Recovering and annotating viral contigs.Viral contigs were recovered from the 6 coassembled biocrust metagenomes using VirSorter 1.0.2, which relies on the detection of “hallmark” viral genes, such as those genes involved in the formation of the viral capsid (72). The Virome db database, without additional viral sequences as references, was used to recover viral contigs from all coassembled biocrust metagenomes. All contigs classified into category 1 or category 2 (complete viral contigs) were considered to be of unambiguous viral origin, while contigs in categories 4 and 5 (proviruses within microbial contigs) were modified to include only the predicted proviral regions, thus excluding host genomic signatures (see host analysis below). Only the viral contigs that scored in the highest confidence categories were retained and analyzed. The viral contigs within our metagenomes are likely sequences derived from double-stranded DNA phages and include prophage elements (temperate phages integrated within microbial chromosomes) and plasmids (extrachromosomal DNA) as well as lytic phages. All predicted viral contigs were retained if they were ≥10 kb in length. The most common class was composed of category 2 viral signatures which were “likely” predictions (n = 1,665), followed by category 1 (“most confident” predictions; n = 184). These viral signals were retrieved from our 28 metagenomes (including the 4 bundle metagenomes), which were coassembled into 5 data sets to assist viral sequence recovery. The removal of short (<10-kb) contigs refined our data set to 1,210 viral genomes or genome fragments. The identified viral protein sequences were clustered and taxonomically classified using vConTACT2 (30), which is supported by the International Committee on Taxonomy of Viruses (ICTV). We built a network on the basis of significant gene-sharing similarities (represented as edges) between viral genomes with contigs represented as nodes. We performed all-vs-all BLASTp queries for comparisons between our 2,140 VirSorted sequences and all predicted viral proteins in the NCBI RefSeq database (v85) with an E value of 1 × 10−3. Protein clusters were defined using the Markov Clustering Algorithm (MCL) and then ClusterONE as described previously (73). These viral genomes represent both high-quality draft genomes (single genome fragment greater than 90% of the expected genome size) and viral genome fragments (single genome fragments less than 90% of the expected genome size or of unknown genome size), defined within the Genomics Standards Consortium framework as the minimum information about an uncultivated virus genome (MIUViG) (16). The UViGs identified by this procedure were visualized within the Cytoscape v3.6.1 program. We also mapped our MAGs and UViGs back to our metagenomes to calculate mapping rates and to determine which genomes were recovered from each time point using Bowtie 2 (74). UViGs were clustered into viral operational taxonomic units (vOTUs) at 95% ANI and 85% AF.
Viral host prediction and CRISPR-Cas analyses.To link viruses with their cognate hosts, we explored the distribution of CRISPR-Cas elements within each annotated MAG and metagenome. Clustered regularly interspaced short palindromic repeats (CRISPR) represent phage infection legacy and thus provide a palpable link between phages and their cognate hosts (36). CRISPRs identified with MinCED 0.3.2 were confirmed using BLASTn queries by allowing for 1 mismatch. We also searched for prophage elements within each MAG using VirSorter, within the CyVerse Discovery Environment.
Auxiliary metabolic gene (AMG) analysis.We considered the 148 viral genomes and genomic fragments for exploration of auxiliary metabolic genes (AMGs) on the basis of VirSorter-predicted contigs (only VirSorter categories 1 and 2, thus excluding proviral contigs). This reduced the possibility of including host genes in our analysis. The identified sporulation related AMGs were reannotated to confirm their classifications through BLASTP comparisons whereby those with E values of 1 × 10−3 and bit scores of >60 were retained. Next, the genomic regions both upstream and downstream of the retained genes were explored for common viral genes, including those encoding terminases, integrases, and viral structural genes. Finally, maximum likelihood trees were built using the AMG sequences and those extracted from the host genomes, as described in the “Phylogenetic tree building” section. The protein structures of our AMGs of interest were predicted with I-TASSER (75), after which we tested for correct predictions using TM-align (76).
Data availability.The sequencing data used in this study are available at the NCBI Sequence Read Archive (project accession number PRJNA395099; sample accession numbers SRX3024638 to SRX3024657; https://www.ncbi.nlm.nih.gov/bioproject/395099).
ACKNOWLEDGMENTS
This work was supported by funds provided by the Office of Science Early Career Research Program, Office of Biological and Environmental Research, of the U.S. Department of Energy under contract no. DE-AC02-05CH11231 to Lawrence Berkeley National Laboratory. Work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under contract no. DE-AC52-07NA27344 (to G.T.).
M.W.V.G. conducted the metagenomic and viromics analyses, formulated ideas, and wrote the manuscript. T.R.N. and T.L.S. conceived the study, and T.R.N., T.L.S., S.R., and G.T. edited the manuscript. T.L.S. collected samples and extracted DNA. S.R. and G.T. assisted with the generation and analysis of viral sequence data. All of us read and approved the final manuscript.
We declare that we have no competing financial interests.
FOOTNOTES
- Received 28 August 2019
- Accepted 7 November 2019
- Published 17 December 2019
This is a work of the U.S. Government and is not subject to copyright protection in the United States. Foreign copyrights may apply.
