Time Series Resolution of the Fish Necrobiome Reveals a Decomposer Succession Involving Toxigenic Bacterial Pathogens

The microbial decomposition of animal tissues is an important ecological process that impacts nutrient cycling in natural environments. We studied the microbial decomposition of a common North American fish (rainbow darters) over four time points, combining 16S rRNA gene and shotgun metagenomic sequence data to obtain both taxonomic and functional perspectives. Our data revealed a strong community succession that was reproduced across different fish and environments. Decomposition time point was the main driver of community composition and functional potential; fish environmental origin (upstream or downstream of a wastewater treatment plant) had a secondary effect. We also identified strains related to the putative pathogen Aeromonas veronii as dominant members of the decomposition community. These bacteria peaked early in decomposition and coincided with the metagenomic abundance of hemolytic toxin genes. Our work reveals a strong decomposer succession in wild-caught fish, providing functional and taxonomic insights into the vertebrate necrobiome.

IMPORTANCE The microbial decomposition of animal tissues is an important ecological process that impacts nutrient cycling in natural environments. We studied the microbial decomposition of a common North American fish (rainbow darters) over four time points, combining 16S rRNA gene and shotgun metagenomic sequence data to obtain both taxonomic and functional perspectives. Our data revealed a strong community succession that was reproduced across different fish and environments. Decomposition time point was the main driver of community composition and functional potential; fish environmental origin (upstream or downstream of a wastewater treatment plant) had a secondary effect. We also identified strains related to the putative pathogen Aeromonas veronii as dominant members of the decomposition community. These bacteria peaked early in decomposition and coincided with the metagenomic abundance of hemolytic toxin genes. Our work reveals a strong decomposer succession in wild-caught fish, providing functional and taxonomic insights into the vertebrate necrobiome. factors affect rivers by introducing pollutants and altering the balance of available nutrients (18)(19)(20), thus influencing fish physiology and their associated microbiomes. Using 16S rRNA sequencing combined with metagenomics, we study both the taxonomic and functional succession of the rainbow darter necrobiome community. We also compare necrobiomes between two different locations in the Grand River (southwestern Ontario, Canada), upstream and downstream of a WWTP, allowing us to analyze community members and their functional potential both spatially and temporally. Studying necrobiome-associated microbial communities provides a unique way to better understand links to aquatic health, fish physiology, and ecosystem dynamics.

RESULTS AND DISCUSSION
Time series community profiling of fish necrobiomes. To examine the structure and temporal succession of aquatic vertebrate necrobiomes, we performed a 16S rRNA-based study of decomposing fish at different time points and locations. We collected female rainbow darters (Etheostoma caeruleum) from the Grand River in Waterloo, Ontario, Canada, both upstream and downstream of the Waterloo wastewater treatment plant (WWTP) (Fig. 1). Individual fish were subjected to decomposition with river water and sediment at room temperature for 1, 4, 8, and 10 days in sterile containers that acted as microcosms of a natural decomposition environment. Sample 16S rRNA gene profiles for fish decomposition microbiomes ("necrobiomes") for these four time points and two water/sediment sources revealed reproducible microbial communities among independent replicates and also between environments (i.e., fish and water source; Fig. 2 and 3). This microbial succession was apparent at the order level of taxonomy (Fig. 2) and at the level of amplicon sequence variants (ASVs) (Fig. 3), although variation in ASV composition was evident among fish samples and environments (Fig. 3). Although infectious lesions began to form on fish sampled on day 1, "bloat-stage" decomposition associated with anaerobic microbial decomposition was not visually apparent until later time points (see Fig. S1b in the supplemental material). The day 1 decomposer communities were composed predominantly of taxa affiliated with Aeromonas and Clostridium, and to a lesser extent by members of the Enterobacteriales. All of these bacterial groups are commonly associated with fish gut microbial communities (21)(22)(23)(24), especially Aeromonas within freshwater fish microbiota (25,26). Species of Clostridium sensu stricto 1 have proteolytic activities and thus may be associated with a more carnivorous diet (24), consistent with the predominantly insectivorous diets of rainbow darters (27). The taxa that were abundant in day 1 fish (e.g., Clostridium sensu stricto 1 sp. ASV 9 and 10, and Aeromonas sp. ASV 3 and 8) decreased in relative abundance over the course of decomposition from 11 to 28% on day 1 to Ͻ 2% on day 10.
Fish sampled on day 4 were associated with anaerobic bloat-stage decomposition and advanced tissue degradation (Fig. S1b). Consistent with previous studies of decomposition using other fish species (14,15), dominant day 4 bacterial phyla detected were Proteobacteria, Firmicutes, and Bacteroidetes (Fig. 2). Compared to day 1 samples, there were considerable changes in degrader community composition for day 4 profiles, with substantial increases in Acetobacteroides ASVs (family Rikenellaceae, order Bacteroidales) from Ͻ 2% on day 1 up to 21% for some ASVs ( Fig. 2 and 3). The Acetobacteroides genus was associated with 45 distinct ASVs across all samples, with 13 ASVs at a relative abundance of Ն2%. Characterized Acetobacteroides members are fermentative, mesophilic, strictly anaerobic, and capable of metabolizing carbohydrates and producing acetate, H 2 , and CO 2 as end products (28). These bacteria also classify under the "Blvii28 wastewater-sludge group" according to the SILVA database. Following day 4, taxa affiliated with Acetobacteroides were the dominant decomposer group, increasing to a relative abundance of as much as 87% in the decomposer community by the final day 10 sampling. Distinct Acetobacteroides ASVs dominated the day 10 decomposer community for fish collected from each of the two river locations (Fig. 3). Acetobacteroides ASV 2 increased to a relative abundance 77% in the fish and water pairing collected upstream of the Waterloo WWTP, whereas Acetobacteroides ASV 1 increased to a relative abundance of 51% in fish collected downstream of the WWTP and left to decompose in water/sediment collected from upstream. These results indicate an influence of the sampled environment on species-or strain-level variations in decomposer communities. The day 4 decomposer communities were also associated with an increase in taxa affiliated with the Selenomonadales order, including Pelosinus and Anaerosinus genera (Fig. 3), which persisted throughout the decomposition process but at low relative abundance. The relative abundance of Selenomonadales increased from 0.0037% on day 1 to 3.8% on day 4, 3.5% on day 8, and 3.3% on day 10.
Influence of spatial location on fish necrobiome succession. Based on sample ordination, necrobiome 16S rRNA gene profiles separated primarily by time point, with distinct microbial communities associated with different stages of tissue decomposition (Fig. 4). Necrobiomes also separated according to sample type, demonstrating distinct profiles for fish that originated upstream versus downstream of the WWTP, but this separation was less apparent than separation based on time. There was no strong effect of water/sediment origin on sample separation. This pattern was found for samples in which fish decomposed in water from the same location and in "swap" experiments in which fish were transferred into sediment/water samples derived from different original locations (e.g., upstream WWTP fish decomposing in downstream WWTP water). Even when isolating the effect of decomposition time, no effect of water/sediment sample origin was detected (Fig. S2). Thus, necrobiomes appear to be influenced primarily by factors occurring prior to decomposition, such as the living fish microbiomes, physiological states, or other fish-environment interactions.
In order to further investigate the influence of fish origin location on decomposer microbial community differences, we calculated differentially abundant taxa (P Ͻ 0.05, Mann-Whitney U test) among fish necrobiomes associated with upstream and downstream WWTP fish sources ( Fig. 5; for additional ASVs, see Table S2 in the supplemental material; for additional statistical methods, see Data Set S1B). Necrobiome-associated , and Pseudomonas (ASV 77). However, these ASVs were low-abundance organisms with a relative abundance of less than 1%. Characterized Peptostreptococcaceae species are anaerobic bacteria that include pathogens associated with tissue infections and antibiotic resistance (29). Known Arcobacter species are aerotolerant and include human and animal pathogens that have been found in groundwater and water reservoirs (30)(31)(32). An additional taxonomic group that increased in relative abundance in necrobiomes originating from fish collected from  downstream of the WWTP was Tolumonas (e.g., ASV 176), a genus with member species associated with the production of toluene from phenol precursors, which are known wastewater pollutants (33,34). Further investigation into the contribution of wastewater effluent on fish decomposition would be helpful to confirm that these differences were indeed linked to WWTP effluent, especially because these putative pathogenic bacteria may contribute to a portion of the energetic costs of fish living downstream of WWTPs (18).
In samples with both fish and water/sediment originating downstream of the WWTP, we also observed higher relative abundance of ASVs associated with Plesiomonas and Lautropia genera (Fig. 3). Again, these ASVs had low (Ͻ1%) relative abundance. Known Plesiomonas species associate with aquatic habitats, cause human infections associated with uncooked shellfish, and have been implicated in infectious outbreaks in regions, including Canada (35,36). Lautropia species have been isolated from the oral cavities of immunocompromised individuals suffering from HIV and cystic fibrosis (37,38).
Metagenomic binning and analysis of decomposition pathways. To explore the genomes and genome-encoded metabolic/functional potential of the necrobiomes, we performed metagenomic sequencing on one replicate for each condition (14 total). Subsequent assembly and binning resulted in four MAGs (metagenome-assembled genomes) with Ͼ85% completion and Ͻ5% redundancy. We examined the taxonomic composition of the MAGs using MetAnnotate (39). These MAGs included two genomes affiliated with Alistipes (Rikenellaceae), a genome annotated as Aeromonas veronii, and a Selenomonadaceae-associated genome ( Table 1). The bins are consistent with ASVs identified by 16S rRNA gene sequencing, corresponding to Acetobacteroides (Rikenellaceae), Aeromonas, and various members of the Selenomonadales (Fig. 2 and 3). Other ASVs identified by 16S rRNA gene sequencing were also recovered in the lower-quality MAGs (Table 1). One bin was affiliated with the genus Pseudomonas, and another bin was affiliated with the family Rikenellaceae.
The relative abundance of Bin_4 (Aeromonas veronii) decreased throughout decomposition from an average relative abundance of 3.7 (day 1) to an average relative abundance of 0.14 (day 10), consistent with our 16S rRNA data (Fig. S3a). Because Aeromonas has been associated with fish gut microbiomes (40)(41)(42)(43)(44), it is possible that Bin_4 and other Aeromonas taxa were initially derived from the fish guts and were important only for early stage decomposition. In contrast, Bin_3 (Rikenellaceae family) may represent a late-stage decomposer because its relative abundance increased in metagenomes from days 8 to 10 of decomposition (average relative abundance of 3.9 on day 8 to an average relative abundance 5.1 on day 10; Fig. S3a). In the downstream fish-upstream sediment/water set, both Rikenellaceae-affiliated bins (Bin_3 and Bin_10) were similar in relative abundance, implying site-specific influences on the relative abundance of different Rikenellaceae-affiliated taxa, consistent with 16S rRNA gene data for Acetobacteroides ASVs (Fig. 3). Phylogenetic analysis of the two Rikenellaceaeassociated bins revealed that Bin_3 was more closely related to Acetobacteroides hydrogenigenes RL-C and Bin_10 was more closely related to Alistipes sp. strain ZOR0009 (Fig. S3b). Bin_9 (Propionispira) was present at low (0.0 to 0.54 average on days 1 to 10; Fig. S3a) relative abundance, close to the sample's mean coverage across the entire course of decomposition, consistent with the abundance patterns seen for Selenomonadales based on 16S rRNA gene data (Fig. 2). Using a KEGG analysis of assembled contigs and binned metagenomes, we examined metabolic pathway potentials associated with decomposition samples. The resulting functional profiles had a highly similar grouping in ordination space compared to the 16S rRNA gene community profiles, whereby samples grouped primarily based on decomposition time point (Fig. S4). Analysis of specific KEGG pathways revealed patterns consistent with a functional succession, mirroring the taxonomic succession described earlier (Fig. 6). Pollutant degradation pathways for polyaromatic hydrocarbons such as naphthalene, styrene, and nitrotoluene showed increased relative abundances on day 1 (13% on average) compared to subsequent time points (6.2% on FIG 6 Selected KEGG pathways displaying significant differential relative abundance across the course of decomposition. Pathways were selected that had an unadjusted P value of Ͻ0.03 after a Kruskall-Wallis test comparing decomposition time (1,4,8 , and 10 days). Shown is the log 10 value of the fractional coverage of the pathway with respect to the total coverage across all the pathways in the sample. Total pathway coverage is also proportionally normalized across every sample. Note that some pathways are based on a few representative genes. For example, coverage of the photosynthesis pathway is mainly derived from genes encoding sodium ion pumps. average). The initial fish bacterial community may have been enriched for microorganisms that could degrade river water contaminants, which can originate from both anthropogenic and natural sources and bioaccumulate in fish (45)(46)(47). Naphthalene degradation in polluted sediment-water systems can be accomplished through several bacterial pathways, and bioremediation of this toxic molecule by native organisms is currently being studied (48)(49)(50). Various biofilm formation pathways were also proportionally abundant (13%) within day 1 metagenomes (Fig. 6), possibly reflecting skin and gut community functions originating prior to decomposition. Degrading river water contaminants and skin and gut biofilm formation may be functions that are more important for the bacterial communities living with their fish host and dealing with possibly contaminated river water than for the necrobiome that formed in our closed system after the fish's death.
Glycan metabolism generally increased in coverage from early stages (2.4% on day 1) to later stages of decomposition (10%). Glycan degradation pathways (e.g., glycosaminoglycans) increased in coverage by days 8 and 10, which may be involved in decomposition of fish skin and intestinal mucins. Late-stage increases in streptomycin, phenylpropanoid, novobiocin, neomycin, kanamycin, and gentamicin biosynthesis pathways (2.4-fold change from day 1 to 10) were also detected, implying that the remaining microorganisms by day 10 possess increased potential for antibiotic synthesis.
These metagenome-wide functional patterns closely matched the functional potentials of individual Aeromonas (early stage) and Rikenellaceae (late stage) bins, when taking into consideration their shifts in relative abundance through the time course (Fig. S5). Genes belonging to pollutant degradation pathways were present in the Aeromonas bin yet mostly absent from other MAGs with lower relative abundance from days 1 and 4 metagenomes. Likewise, biofilm formation pathway genes had a 6.2-foldhigher frequency in the Aeromonas bin compared to the Acetobacteroides/Alistipes bins. In contrast, antibiotic biosynthesis pathway genes had a 2.5-fold-higher frequency in the Rikenellaceae-associated bins, in addition to multiple key glycan degradation genes. Thus, the detected shifts in functional profiles were in part due to the hand-off microbial community dominance from Aeromonas to Rikenellaceae. It is important to note that these apparent late-stage functional shifts could also be important for earlier phases when Rikenellaceae initially began to increase in relative abundance.
Our data suggest strong Acetobacteroides dominance in late-stage rainbow darter necrobiomes ( Fig. 2 and 3). Because related species have been implicated in anaerobic sugar fermentation (28), we investigated the two MAGs affiliated with these bacteria for glycolytic enzymes. Both Bin_3 and Bin_10 possess a complete glycolysis pathway as well as L-lactate dehydrogenase for anaerobic fermentation (Fig. S6). Bin 3 genes also encode pyruvate dehydrogenase, aldehyde dehydrogenase, and enzymes for conversion of D-fructose, D-fructose-1-phosphate (D-fructose-1P), and D-mannose-6P to glycolysis precursors. Based on a previous analysis of decomposition pathways (51), Bin_3 and Bin_10 genes also encode components of potential pathways for production of indole (EC 4. Previous research showed that Acetobacteroides hydrogenigenes RL-C can produce acetate and carbon dioxide from glucose fermentation (28). Because these metabolites could potentially be converted to methane by methanogenic archaea (52), we analyzed the assembled metagenomic data for archaea-associated contigs. Contigs taxonomically affiliated with methanogenic archaea were identified from almost all samples, including species of the classes Methanobacteria, Methanococci, and Methanomicrobia (Fig. S7a). A 1.7-fold-higher relative abundance of these contigs was observed after day 1 (Fig. S7a), indicating that methanogenic archaea may have increased in relative abundance early in decomposition, coinciding with anoxic conditions and the generation of acetate or carbon dioxide by Rikenellaceae bacteria. However, no archaeal taxa were identified in the 16S rRNA. This potentially reflects an increased ability to detect Time Series Metagenomic Analysis of a Fish Necrobiome low-abundance archaeal organisms in our metagenomic data set, perhaps due to the availability of more taxonomic markers and/or lower archaeal 16S copy numbers (53).
A toxigenic strain of Aeromonas veronii is a dominant member of the necrobiome. Because Bin_4 affiliated with A. veronii, a well-established pathogen of fish and humans (36,(54)(55)(56)(57)(58)(59)(60), and a common inhabitant of the fish gut microbiome (40)(41)(42)(43), we explored its phylogenetic position, functional profile, and virulence repertoire. A maximum likelihood phylogeny of A. veronii and other related Aeromonas genomes from the NCBI was constructed based on a concatenated alignment of conserved ribosomal marker genes (Fig. 7a). Within this phylogeny, Bin_4 grouped with a clade of A. veronii genomes but as a basal lineage outgrouping all A. veronii species except AMC34.
We used the VFanalyzer from the Virulence Factor Database (VFDB) to detect virulence factors within Bin_4 and compared it to a reference Aeromonas strain, A. veronii B565. This bin contained virulence-related genes for adherence, iron uptake, and secretion systems (Data Set S1C). Indeed, a total of 54 genes that were associated with secretion systems were identified, compared to only 15 in A. veronii B565. In addition, we identified 13 genes associated with endotoxin production. Like A. veronii B565, Bin_4 genes encoded hemolysin III, hemolysin HlyA, and a thermostable hemolysin gene (Fig. 7b). We also recovered a relatively small incomplete bin (Bin_11, 0.64 Mb, 717 coding sequences [CDSs], 321 contigs) that correlated with Bin_4 in relative abundance (Fig. S7b). This small bin affiliated with Aeromonas veronii and also included a gene encoding aerolysin toxin production (Data Set S1D). Based on metagenomic mapped read coverage, the relative abundance of genes encoding Aeromonas toxins increased on day 4 of decomposition (Fig. 7c), indicating an enrichment in Aeromonas strains carrying hemolytic proteins. A possible explanation for this is that lytic toxins, including those from Aeromonas, may function in host cell lysis during decomposition and therefore peak in relative abundance during earlier stages of decomposition. Bin_4 also possessed genomic potential for decomposition-related pathways, including histidine degradation (contains EC 4. Conclusion. Overall, our microcosm study of rainbow darter fish decomposition revealed a highly reproducible microbial succession throughout the time course, even across different fish and water/sediment sources. The location of the fish when sampled (upstream or downstream of the WWTP) also affected its decomposition profile, suggesting that necrobiomes may be influenced by prior fish-environment interactions. Together, our data suggest that environmental interactions may shape the initial gut community and/or the physiological state of the fish, which then seeds or impacts the later necrobiome community and its succession.
Both 16S and metagenomic analysis revealed a strong succession in which initial time points were dominated by Clostridiaceae and Aeromonas, with Rikenellaceae species appearing by day 4 and becoming major community members by day 10. Analysis of functional profiles inferred from the metagenomic data revealed common decomposition pathways, as well as temporal shifts in function that mirrored taxonomic succession. Notably, pollutant degradation pathways and biofilm formation pathways were enriched in the early stages of decomposition and associated with Clostridiaceae and Aeromonas, and glycan metabolism and antibiotic synthesis increased in later stages and associated with Rikenellaceae. Last, we identified a toxigenic Aeromonas strain that was a dominant member of the necrobiome community. The presence of numerous hemolytic toxin genes in this organism suggests a potential role for toxins in the decomposition of host tissues as proposed previously (11). Further work investigating the prevalence and function of toxigenic bacterial species in decomposer communities will be important to explore their broader ecological roles and niches within natural ecosystems.

MATERIALS AND METHODS
Fish collection. On 24 October 2016, female rainbow darters (Etheostoma caeruleum) were collected from the Grand River (Fig. 1 (43°29=16==N; 80°30=25==W). Forty-two fish (21 from each site) were collected using a backpack electrofisher (Smith Root, LR-20) and euthanized quickly with a sharp blow to the head. Then each fish was placed in an autoclaved 250-ml mason jar microcosm that contained a mixture of water and river substrate (see Table S1 for river water quality metadata and Fig. S1a for an example mason jar setup). The lids were closed, but not sealed, in order to ensure oxic conditions that would accompany natural in-river decay events. The jars were then left to decay in a fume hood at room temperature. Three samples containing both fish and water/sediment from the same site were left to decompose for 1 day (24 h), 4 days, 8 days, and 10 days for both the WMR and EIT sites, totaling 24 fish. For additional treatments to assess differences in water quality and aquatic microorganisms, three samples containing fish and water/sediment from different sites (i.e., WMR fish in EIT conditions and EIT fish in WMR conditions) were allowed to decay for 4, 8, and 10 days, totaling 18 fish. At each time point, decay was documented (Fig. S1b), and fish were removed from the replicate jars, then rinsed with sterile water, and ground with liquid nitrogen using a clean mortar and pestle. The powdered tissue was stored at -80°C prior to genomic DNA extraction. Following this lysis incubation, 700 l of the lysate was extracted with an equal volume of phenol and centrifuged at 10,000 ϫ g for 5 min. The aqueous phase was retained and twice extracted with equal volumes of phenol-chloroform-isoamyl alcohol (25:24:1), followed each time with centrifugation at 10,000 ϫ g for 5 min. One volume of isopropanol was used to precipitate aqueous phase DNA in a new ultracentrifuge tube, followed by centrifugation at 13,000 ϫ g for 10 min at room temperature. The resulting pellet was washed twice with 70% ethanol, dried, and then dissolved in 50 l of DNase-and RNase-free H 2 O (Sigma) at 50°C for 15 min. The quantity and quality of DNA were determined with a SpectraDrop (Molecular Devices) and stored at -20°C prior to sequencing.
16S rRNA gene and metagenomic sequencing. Extracted DNA was amplified in triplicate using Pro341F and Pro805R universal prokaryotic primers (61). Triplicate amplicons were pooled, gel quantified, and sequenced to a depth of at least 30,000 paired-end reads per sample using the MiSeq reagent kit v3 (2 ϫ 300 cycles; Illumina).
For metagenomic sequencing, genomic DNA (1 ng) was fragmented and individually barcoded using the Nextera XT DNA Library Prep kit (Illumina) following the supplier's guidelines. Small fragments of library DNA were removed by adding 0.6 volumes of AMPure XP beads (Beckman Coulter). After washing twice with 80% ethanol and air drying for 10 min, DNA was eluted from the beads with 10 mM Tris-HCl (pH 8.5). Purified library DNA was quantified with the Qubit dsDNA (double-stranded DNA) HS (highsensitivity) assay kit, diluted to 4 nM with the Tris-HCl buffer and then pooled in an equal volume. Library DNA was denatured with equal volumes of 0.2 N NaOH, diluted to 7 pM with hybridization buffer HT1, and sequenced with MiSeq reagent kit v2 (2 ϫ 250 cycles; Illumina).
16S rRNA gene analysis. Demultiplexed sequences were processed using DADA2 v1.4 (62), managed through QIIME2 v.2017.10 (63). Briefly, forward and reverse reads were truncated with decreasing quality metrics while maintaining sequence overlap (ϳ250 bases). Primers were removed, and paired reads were assembled after error modeling and correction, creating amplicon sequence variants (ASVs). Chimeric ASVs were removed by reconstruction against more abundant parent ASVs. The resulting ASV table was constructed for downstream analysis (see Data Set S1A in the supplemental material).
Taxonomy was assigned to representative sequence variants using a naive Bayesian classifier implemented in QIIME2 with scikit-learn (v.0.19.0), trained against SILVA release 128 (64), clustered at 99% identity, and trimmed to the amplified region. Assignments were accepted above a 0.7 confidence threshold.
For ordination, we used a proportion matrix of ASVs across each sample with a sparsity cutoff (i.e., ASV detected in at least 3 of 42 samples). The metaMDS() and envfit() scripts from vegan package v2.4-2 in R were used to calculate ordination coordinates and data vectors. A stress or Shepard diagram was generated with stressplot() from the vegan package to determine the nonmetric fit. The ASVs with significant rank sum differences in sample proportion were calculated with the Mann-Whitney test in R. Multiple hypothesis correction of P values was performed using the p.adjust() function in R with the Benjamini-Hochberg model. We also calculated differential taxon relative abundance using a variety of methods (metagenomeSeq, edgeR, DESeq2, and LEfSe) as implemented in the Marker Data Profiling pipeline from MicrobiomeAnalyst (65) with default settings on 20 March 2020.
For metagenomic and bin functional analysis, KEGG (Kyoto Encyclopedia of Genes and Genomes) annotations were identified with GhostKOALA (67). The average coverage for each gene (per base pair), normalized by dividing by the average sample coverage (per base pair), was summed to give a total coverage value for each KEGG pathway. The decostand() function from the vegan package v2.4-2 in R was used to determine the fractional value of each pathway with respect to the total summed coverage across all KEGG pathways detected in the sample. A Kruskall-Wallis test was done in R to identify KEGG pathways with significantly different distributions by day of decomposition. The decostand() function was also used to proportionally normalize each pathway value across every sample for plotting. For the bin functional analysis, the frequency of each KEGG orthology (KO) annotation in each MAG bin was counted. These counts were summed for each KEGG pathway, and fractional values were calculated across all KEGG pathways detected in the bins as before.
The VFanalyzer software from the Virulence Factor Database (VFDB) (68) identified virulence factors in the predicted coding sequences of Bin_4 using Aeromonas veronii B565 as a representative genome. We also used the domain architecture from the Aeromonas toxin gene set from the VFDB to identify Aeromonas toxin genes in the coassembly. Putative toxins longer than 150 amino acids were assessed with BLASTp for Aeromonas taxonomy and gene annotation. Data availability. All 16S rRNA gene and metagenomic sequencing data for this project were deposited into the NCBI Short Read Archive (SRA) under BioProject accession no. PRJNA604775.

SUPPLEMENTAL MATERIAL
Supplemental material is available online only.