Temporal dynamics in microbial soil communities at anthrax carcass sites

Abstract

Background: Anthrax is a globally distributed disease affecting primarily herbivorous mammals. It is caused by the soil-dwelling and spore-forming bacterium Bacillus anthracis. The dormant B. anthracis spores become vegetative after ingestion by grazing mammals. After killing the host, B. anthracis cells return to the soil where they sporulate, completing the lifecycle of the bacterium. Here we present the first study describing temporal microbial soil community changes in Etosha National Park, Namibia, after decomposition of two plains zebra (Equus quagga) anthrax carcasses. To circumvent state-associated-challenges (i.e. vegetative cells/spores) we monitored B. anthracis throughout the period using cultivation, qPCR and shotgun metagenomic sequencing.

Results: The combined results suggest that abundance estimation of spore-forming bacteria in their natural habitat by DNA-based approaches alone is insufficient due to poor recovery of DNA from spores. However, our combined approached allowed us to follow B. anthracis population dynamics (vegetative cells and spores) in the soil, along with closely related organisms from the B. cereus group, despite their high sequence similarity. Vegetative B. anthracis abundance peaked early in the time-series and then dropped when cells either sporulated or died. The time-series revealed that after carcass deposition, the typical semi-arid soil community (e.g. Frankiales and Rhizobiales species) becomes temporarily dominated by the orders Bacillales and Pseudomonadales, known to contain plant growth-promoting species.

Conclusion: Our work indicates that complementing DNA based approaches with cultivation may give a more complete picture of the ecology of spore forming pathogens. Furthermore, the results suggests that the increased vegetation biomass production found at carcass sites is due to both added nutrients and the proliferation of microbial taxa that can be beneficial for plant growth. Thus, future B. anthracis transmission events at carcass sites may be indirectly facilitated by the recruitment of plant-beneficial bacteria.

Keywords: Bacillus anthracis; Metabolism; Metagenomics; Microbial diversity; Semi-arid; Shotgun sequencing; Sporulation; Taphonomy; Time-series analysis.

Fig. 1

Estimation of B. anthracis abundance in soil samples using qPCR, metagenomic reads mapping and cultivation. For all panels, sampling time-points are on the x-axis. a and d qPCR results with estimated number of B. anthracis genomes per gram soil along the y-axis at different time-points (x-axis) for Carcass 1 (Ca1) (a) and Carcass 2 (Ca2) (d). The spiked sample is control soil from day 0 with 3.4 × 10⁶ Stern34F2 B. anthracis spores added. P1 and P2 represent the two replicates at each time-point. b, c, e and f shows the results from mapping the metagenomic reads against the B. anthracis isolate genomes (K1 /K2) isolated from Ca1 and Ca2 as well as the other reference strains; B. cereus E33L, B. thuringiensis HD-771 and B. subtilis 168. Reads matching the references are on the y-axis. b and e and (c and f) shows the mapping of the metagenomes against the K1 and K2 strain, respectively as well as the other reference strains. Ca1 in panes (B and C) and Ca2 in panes (E and F). (G) B. anthracis spore counts of soil samples from Ca1 and Ca2 in spores per gram dry soil. Spore counts were estimated by culturing of heat-shocked soil samples and adjusted to dry weight of soil

Fig. 2

OTU richness of 16S /18S rRNA sequences from soil metagenomic samples. 16S /18S rRNA sequences were extracted with metaxa2 and diversity was analysed with MetaAmp. a Rarefaction curves. X-axis indicates number of rRNA sequences, and y-axis shows OTU numbers b Rank abundance of OTUs, where OTU rank is on the x-axis and the abundance per OTU is on the y-axis. The solid lines are Ca1 samples and the dotted lines are Ca2 samples. Line colour indicates time-points as indicated in the legend

Fig. 3

Soil microbial community relationships by time-point and carcass site. Principal coordinate analysis of a distance matrix created by normalised total counts and using Bray-Curtis dissimilarity. Ca1 is represented by the circles and Ca2 by the triangles, the time-points are visualised in different colours and the arrows are pointing in the direction of increasing days, red arrows for Ca1 and light blue for Ca2. The dotted arrows show the relationship between the day 30 to the Ctrl0 sample

Fig. 4

Temporal dynamics of microbial order abundances at carcass 1 and 2. a Heatmap visualisation of the relative abundance for the 50 most abundant orders at each time-point for carcass 1 and 2. b Barplots that show the log5 average fold change of the raw reads counts across the time-series for the 50 taxa shown in figure (a). The taxa in the heatmaps (a) are sorted using the order abundances in the control sample of Carcass1 at day 0 (Control). In order to visualise the relative abundances (a) values were log5 normalised and then scaled so that the sum of each column equals 1. Eukaryotic orders are marked with *, the remaining orders are bacterial. The bottom seven entries have not been classified to lower phylogeny than Kingdom, Phylum or Class

Fig. 5

Heatmap of significantly different metabolic pathways per time-point. KEGG metabolic pathways significantly correlating with AGS for microbial communities of Ca1 and Ca2. Pathways were determined by MEGAN classification. KEGG pathways abundances were normalised with DESEQ2 [96] and correlated to the AGS using the Spearman correlation method with the False Discovery Rate (FDR) test to calculate probabilities (FDR cut-off at 0.05). KEGG-pathways correlating positively or negatively with a p-value <0.05 are shown for Ca1 (blue) and Ca2 (red) for all the time-points, the AGS is shown in the top panel. # indicates pathways that are significant in both Ca1 and Ca2, § indicates pathways that are significant in Ca2, and the rest are only significant in Ca1. Pathway abundances were centred and scaled per row and positive and negative correlations were clustered based on the sign of the Spearman correlations