Playing with FiRE: A genome resolved view of the soil microbiome responses to high severity forest wildfire

Warming climate has increased the frequency and size of high severity wildfires in the western United States, with deleterious impacts on forest ecosystem resilience. Although forest soil microbiomes provide a myriad of ecosystem functions, little is known regarding the impact of high severity fire on microbially-mediated processes. Here, we characterized functional shifts in the soil microbiome (bacterial, fungal, and viral) across wildfire burn severity gradients one year post-fire in coniferous forests (Colorado and Wyoming, USA). We generated the Fire Responding Ecogenomic database (FiRE-db), consisting of 637 metagenome-assembled bacterial genomes, 2490 viral populations, and 2 fungal genomes complemented by 12 metatranscriptomes from soils affected by low and high-severity, and complementary marker gene sequencing and metabolomics data. Actinobacteria dominated the fraction of enriched and active taxa across burned soils. Taxa within surficial soils impacted by high severity wildfire exhibited traits including heat resistance, sporulation and fast growth that enhanced post-fire survival. Carbon cycling within this system was predicted to be influenced by microbial processing of pyrogenic compounds and turnover of dominant bacterial community members by abundant viruses. These genome-resolved analyses across trophic levels reveal the complexity of post-fire soil microbiome activity and offer opportunities for restoration strategies that specifically target these communities.


Introduction
additionally pre-treated with solid-phase extractions using Agilent Bond Elut-PPL 3 mL columns 118 and diluted to 50 ppm (Agilent Technologies, DE, USA) following standard lab protocol 29 . A 12 119 Tesla (12T) Bruker SolariX FTICR-MS located at the Environmental Molecular Sciences 120 Laboratory in Richland, WA, USA was used to collect DOM high-resolution mass spectra from 121 each DOM sample. Samples were directly injected into the instrument using a custom automated 122 direction infusion cart that performed two offline blanks between each sample and using an Apollo 123 II electrospray ionization (ESI) source in negative ion mode with an applied voltage of -4.2kV. 124 Ion accumulation time was optimized between 50 and 80 ms. One hundred and forty-four 125 transients were co-added into a 4MWord time domain (transient length of 1.1 s) with a spectral 126 mass window of m/z 100-900, yielding a resolution at m/z 400. Spectra were internally 127 recalibrated in the mass domain using homologous series separated by 14 Da (CH2 groups). The 128 mass measurement accuracy was typically within 1 ppm for singly charged ions across a broad 129 m/z range (100 m/z -900 m/z). Bruker Daltonics DataAnalysis (version 4.2) was used to convert 130 mass spectra to a list of m/z values by applying the FTMS peak picking module with a signal-to-131 noise ratio (S/N) threshold set to 7 and absolute intensity threshold to the default value of 100. 132 Chemical formulae were assigned with Formularity 30 based on mass measurement error < 0.5 ppm, 133 taking into consideration the presence of C, H, O, N, S and P and excluding other elements. This 134 in-house software was also used to align peaks with a 0.5 ppm threshold. The R package 135 ftmsRanalysis 31,32 was then used to remove peaks that either were outside the desired m/z range 136 (200 m/z -900 m/z) or had an isotopic signature, calculate nominal oxidation state of carbon 137 (NOSC), and assign Van Krevelen compound classes. Raw FTICR-MS data is provided in archive 138 (doi:10.5281/zenodo.5182305).

140
Kendrick mass defect (KMD) analysis and plots were employed to identify potential increasing 141 polyaromaticity across the burn severity gradient. The KMD analysis was done using the C4H2 142 base unit (50 amu) to represent the addition of benzene to a separate molecular benzene. Research, CA, USA). Ribosomal RNA was removed from total RNA and libraries were prepared 234 using the Takara SMARTer Stranded Total RNA-Seq Kit v2 (Takara Bio Inc, Shiga, Japan).

235
Samples were sequenced on the NovaSEQ6000 platform on a S4 flow cell using 151-bp paired-236 end reads at Genomics Shared Resource, Colorado Cancer Center, Denver, CO, USA. Adapter 237 sequences were removed from raw reads using Bbduk (https://jgi.doe.gov/data-and-238 tools/bbtools/bb-tools-user-guide/bbduk-guide/) and sequences were trimmed with Sickle v1.33 42 . 239 Trimmed reads were mapped to metagenome assemblies using BBMap (parameters: 240 ambiguous=random, idfilter=0.95; v38.70). Mappings were filtered to 95% identity and counts 241 were generated using HTSeq 55 . For broad analysis of differential expression, the dataset was 242 filtered to only transcripts which were successfully annotated by DRAM (n=146,895) and 243 DESeq2 56 was used to identify transcripts that were differentially expressed in either burn severity 244 within soil horizons (e.g., high vs. low in organic horizon soils and vice versa). The same analysis 245 was also run on the combined HMM output described above (1,189 total transcripts). We 246 normalized our dataset by calculating the gene length corrected trimmed mean of M values 57 247 (geTMM) using edgeR 58 to normalize for library depth and gene length. To identify transcripts 248 that were highly expressed in any given condition, we filtered the data to only transcripts that were 249 in the upper 20% of TMM for 2 of the 3 samples in any one condition ( prior studies 17,20,65,66 , bacterial communities in burned soils were characterized by lower diversity 295 and were enriched in Actinobacteria and Bacteroidetes relative to unburned controls (Figure 1). 296 Specifically, the Actinobacteria genera Arthrobacter, Modestobacter, Blastococcus, and 297 Actinomadura had the largest relative abundance increases in burned soils relative to control soils 298 and were discriminant features of these conditions (Supplementary data 8). The diversity of 299 fungal communities in surficial soils also decreased with fire ( Figure 2E) and, similar to previous 300 studies 18,67,68 , shifted from Basidiomycete-to Ascomycete-dominated with the Basidiomycota 301 relative abundance decreasing by ~58% (Figure 1). Discriminant fungal taxa included ASVs from 302 the Sordariomycetes, Saccharomycetes, and Dothideomycetes, taxa also found in previous fire 303 studies 66,68 (Supplementary data 8).

305
While shifts in community composition between burned and control conditions were observed in 306 both organic and mineral soil horizons, the surficial organic layers changed more than the deeper 307 mineral horizon soil. Microbial diversity generally decreased with increasing fire severity in the 308 organic horizons, although differences between moderate and high severity were statistically 309 indistinct (Figure 2). Similarly, as fungal and bacterial diversity decreased with burn severity, beta 310 dispersion ('distance to centroid') calculations revealed increasingly similar microbial community 311 structures ( Figure S3) with a less complex bacterial and fungal community structure (Table S2). 312 These shifts resulted in significant dissimilarity between microbial communities in organic 313 horizons impacted by either low or high severity wildfire (bacterial ANOSIM R=0.15, p<0.05; 314 fungal ANOSIM R=0.25, p<0.05). In contrast, mineral soils that were less impacted by wildfire 315 displayed an opposite effect, with increasing beta dispersion after wildfire that signifies greater 316 bacterial community dissimilarity ( Figure S3). Stochastic community shifts in deeper soils may 317 follow wildfire, potentially due to spatially heterogeneous changes in soil chemistry and nutrient 318 availability. Together, these data highlight the susceptibility of highly combustible surficial 319 organic soil horizons to wildfire, resulting in less diverse and inter-connected microbial 320 communities. In contrast, mineral soils are likely more insulated from the effects of fire, and 321 therefore the soil microbiome displays a more muted response to wildfire effects .  322  323  324  325  326  327  328  329  330  331  332  333  334  335  336  337  338  339  340  341  342  343 344 345 Figure 1 The percent change in relative abundance from control to low, moderate, and high 346 severity in organic soil horizons of each main bacterial and fungal phylum. Phyla with relative 347 abundance less than 0.5% were discarded for this analysis. Note that although the Firmicutes have 348 the largest increase with burn (inset) their overall relative abundance is still low.  High severity wildfire exerts a pulse disturbance on surficial organic soil horizons by exposing 396 them to high heat. Thus, the taxa that constitute the post-fire microbiome and fill new niches in 397 the soil likely encode traits for thermal resistance. Nearly all the aforementioned Actinobacterial 398 MAGs encoded sporulation genes, indicating that spore formation is likely a key strategy to ensure 399 survival and colonization post-fire ( Figure 3A). These MAGs also all contained genes encoding 400 heat shock proteins and molecular chaperones which may further facilitate thermal resistance 401 (Figure 3A). In the majority of these MAGs, thermal resistance was complimented by genes for 402 mycothiol biosysnthesis, a compound produced by Actinobacteria that aids in oxidative stress 403 tolerance 72 . Finally, MAG RYN_93 encoded ectB for synthesizing ectoine, a compatible solute for 404 environmental stress tolerance 73 . In general, MAGs recovered from High O samples also had 405 significantly higher GC content than MAGs from mineral soil layers (e.g., A horizons) impacted 406 by low ('Low A') and high ('High A') wildfire severity ( Figure S5) ( Figure 3B). These insights suggest that abundant bacteria sampled one year post-wildfire likely 426 occupied niches in the immediate aftermath of wildfire through rapid cell growth. In contrast, these 427 patterns were absent from MAGs from High A soils ( Figure S7)  Deeper burned soils are dominated by distinct microbial membership 462 High A soils hosted a greater diversity of enriched and active MAGs relative to the overlying 463 organic soil horizons. Actinobacteria again contributed strongly to these signals, with the 464 dominant, active MAGs in the mineral horizons affiliated with the families Streptosporangiaceae, 465 Solirubrobacteraceae, Frankiaceae, and Streptomycetaceae (MAGs RYN_173, RYN_225, 466 RYN_228, RYN_220, RYN_230, RYN_265). Additional highly abundant and active MAGs in 467 High A samples were affiliated with the Eremiobacterota, Dormibacterota, and Proteobacteria 468 phyla (MAGs RYN_132, RYN_309, RYN_342, RYN_347, RYN_607). Together, these MAGs 469 accounted for ~20% of the High A community composition, and nearly 14% of total MAG gene 470 expression. These MAGs were also active in Low A soils, albeit to a lesser extent (accounting for 471 ~5% of total expression).

473
Genes associated with thermal resistance were again common in the dominant High A MAGs.

474
Three of the MAGs discussed here (RYN_220, RYN_225, RYN_342) had metabolic potential to 475 synthesize the environmental stress protectant mycothiol 72 . The second most active MAG in High 476 A samples (RYN_220) was affiliated with the Streptomyces which have been shown to form 477 sporogenic structures (aerial hyphae) in response to adverse conditions (i.e., high temperatures, 478 nutrient depletion) 81 . We observed an overall increase in spore-forming gene coverage from Low 479 A to High A conditions (Figure 2A) and all MAGs discussed here encoded at least one sporulation 480 gene. Additionally, members of Streptomyces are known for their ability to scavenge nutrients and 481 utilize a diverse array of organic substrates 82 , which is evidenced here through expression of 482 abundant and diverse CAZymes that are predicted to target chitin, polyphenolics, and pectin, 483 among other substrates ( Figure S6).

485
Associated with the enrichment and activity of Streptomyces in High A soils, we observed high 486 expression of genes encoding the biosynthesis of cobalamin (vitamin B12, cob genes). Cobalamin 487 production is conserved within a relatively small group of microorganisms -including 488 Streptomyces -and can serve as a keystone function within ecosystems 83 . The entire aerobic 489 cobalamin biosynthesis pathway was expressed in Low A and High A samples (Figure S8; 490 pathway adapted 83,84  wildfire enriching taxa that encode this trait (i.e., Streptomyces). Given the noted importance of 500 this cofactor in mediating a range of critical soil microbiome functions, this process could 501 potentially aid in plant reestablishment 85 , enhance ecosystem function across trophic levels 86 , and 502 thus be beneficial for the overall restoration of post-fire landscapes. 503 504 Actinobacteria catalyze degradation of pyrogenic organic matter 505 During wildfire, SOM may be transformed to increasingly aromatic molecular structures that are 506 commonly considered less labile for microbial utilization 87 . Mass spectrometry DOM analyses 507 found evidence for increasing aromaticity with burn severity in organic horizon soils (Figure 4). 508 Low severity wildfire drives an increase in DOM aromaticity ( Figure 4B) but also an 509 accumulation of other unique compounds likely resulting from incomplete combustion of SOM 88 510 ( Figure 4A). In contrast, the unique compounds formed following moderate and high severity 511 wildfire were constrained to the aromatic region ( Figure 4A). Reflecting the more insulated, 512 deeper mineral soils, these aromaticity index trends were not identified in DOM released from 513 mineral horizons ( Figure S9). To estimate lability of these DOM pools, we calculated NOSC 514 which can reveal the potential thermodynamic favorability of a carbon substrate, with higher 515 NOSC values theoretically yielding a lower ΔG C ox (i.e., more favorable) when coupled to the 516 reduction of an electron acceptor 89 . Unique formulas detected in High O samples had significantly 517 higher NOSC values than both Low O and control samples, indicating increasing thermodynamic 518 favorability of DOM following severe wildfire ( Figure 4C). 519 520 Wildfire severity-driven chemical changes likely mediate the extent of microbial DOM processing. 521 We focused on the ability of community members to process catechol and protocatechuate, the 522 two intermediate products formed during aerobic degradation of diverse aromatic compounds 90 . 523 The genomic potential for these reactions was broadly distributed across burn severities and was 524 dominated by Actinobacteria and Proteobacteria ( Figure 5); 80 and 226 MAGs encoded at least 525 50% of the catechol and protocatechuate ortho-cleavage pathways, respectively, including most of 526 the featured High O and High A MAGs ( Figure 5C). Meta-cleavage pathways were also broadly 527 represented within the MAG database ( Figure S10). There was additional evidence for activity of 528 both pathways regardless of soil horizon and wildfire severity ( Figure 5A). In High O samples, 529 the Arthrobacter MAG RYN_101 alone was responsible for ~44% of catA (catechol 1,2-530 dioxygenase) gene expression, and therefore likely plays a key role in catechol degradation. 531 Contrastingly, in High A samples, the Streptosporangiaceae MAG RYN_225 was responsible for 532 ~46% and 23% of the expression of pcaGH (protocatechuate 3,4-dioxygenase) and pcaC (4-533 carboxymuconolactone decarboxylase) that catalyzes protocatechuate degradation ( Figure 5C). 534 However, none of the dominant MAGs across High O and High A samples encoded the entire 535 catechol or protocatechuate ortho-cleavage pathway (Figure 5C), indicating that metabolic 536 handoffs between community members are likely important for complete degradation of these 537 compounds. 538  545  546  547  548  549  550  551  552  553  554   555  556  557  558  559  560  561  562  563  564  565  566  567  568  569  570  571  572  573  574  575  576  577  578  579  580  581  582  indicating the log normalized sum TMM of the gene for high severity organic and mineral soil 588 horizons. Asterisks indicate genes that are differentially expressed in the condition (p<0.05). (C) 589 The genomic potential and expression of each gene in the pathway for the MAGs enriched and 590 active in high severity organic and mineral soil samples. The above bar chart shows the featured 591 MAG relative abundance in that condition, colored by featured condition. 592 593 594 Viral dynamics impact community structure and function of the burned soil microbiome 595 We recovered 2399 distinct DNA viral populations and 91 distinct RNA viral populations 596 (vMAGs) from the metagenomic and metatranscriptomic assemblies. Highlighting the viral 597 novelty in these ecosystems, 945 of the DNA and RNA vMAGs were previously undescribed (only 598 clustering with other vMAGs from this study) and 92 were given taxonomic assignments, with the 599 majority (n=86) within the Caudovirales order (Table S5). Both DNA and RNA viral communities 600 mirrored beta diversity trends observed in the bacterial and fungal communities; those in mineral 601 soil horizons were less constrained in multivariate space compared to communities in organic 602 horizons, further highlighting the homogenizing influence of wildfire ( Figure S11). Additionally, 603 DNA and RNA viral community composition was not significantly different between low-and 604 high-severity impacted soils (ANOSIM R = 0.007 and -0.12, respectively; p > 0.1) but significantly 605 differed between the two soil horizons (ANOSIM R = 0.59 and 0.57, respectively; p<0.05).

607
Given the importance of viral activity on soil microbiome community composition and 608 function 91,92 , we identified potential virus-host linkages that could offer insights into how viruses 609 target bacteria. Many abundant and active MAGs (n=94) -including 32 from the Actinobacteria 610 -encoded CRISPR-Cas arrays with an average of ~18 spacers (max 210 spacers; Supplementary  611 data 3). Through the matching of protospacers to sequences in vMAGs, we linked 9 vMAGs with 612 4 bacterial hosts (RYN_115, RYN_242, RYN_436, RYN_542) from the phyla Actinobacteria, 613 Planctomycetota, and Proteobacteria. While each of these four MAGs were active (via transcript 614 expression), the RYN_242 MAG (Solirubrobacteraceae) was among the top 3% of active MAGs 615 across all four conditions, suggesting that viruses are targeting active bacteria. We expanded upon 616 potential virus-host linkages using VirHostMatcher 64 (d2 * value < 0.25), revealing higher numbers 617 of viral linkages with more abundant host MAGs (Figure 6). For example, the MAGs of interest 618 from High O and High A samples had above average numbers of putative viral linkages (average 619 of 196; Figure 6). Moreover, there were 198 vMAGs with putative linkages to all 8 featured 620 Actinobacteria MAGs from High O soils, potentially due to conserved nucleotide frequencies. The 621 shared 198 vMAGs made up ~11% of the viral community in High O samples, again suggesting 622 that the most abundant and active bacteria in burned soils are being actively targeted by abundant 623 phage, potentially impacting soil C cycling via release of labile cellular components following cell 624 lysis. There is also evidence in this system for the "piggyback-the-winner" viral strategy, where 625 lysogenic lifestyles are favored at high microbial abundances and growth rates 93 . Of our 2399 626 DNA vMAGs, we identified 185 with putative lysogenic lifestyles based on gene annotations for 627 integrase, recombinase, or excisionase genes, 32 of which had nucleotide frequency-based 628 linkages to all the featured Actinobacteria MAGs in High O samples.

630
To further investigate potential viral roles in C cycling in burned soils, we characterized the 631 putative AMG repertoire of the vMAGs. Viruses utilize AMGs to "hijack" and manipulate the 632 metabolisms of their host; one study in thawing permafrost soils found AMGs associated with 633 SOM degradation and central C metabolism, suggesting that viruses can play a direct role in 634 augmenting the soil C cycle 92 . There were 773 total putative AMGs detected in 445 vMAGs within