Postmortem samples included 28 male and 17 female corpses from the Alabama Department of Forensic Sciences in Montgomery, AL and The Office of the District One Medical Examiner in Pensacola, FL, United States. Demographic data were collected on each of the 45 corpses (i.e., age, gender, ethnicity, cause of death, PMI) (Supplementary Table S1). The study was approved by Alabama State University Institutional Review Board (IRB) number 2016011. Time of death of each corpse was certified from official Daily Crime Logs created by local police departments. Bodies were stored in the morgues at 1°C until time of tissue dissection. Approximately 10 mg of liver and spleen tissues were dissected using sterile scalpels and placed in polyethylene bags in an examination area at 20°C. Organs were transported on dry ice to the Thanatos Laboratory at Alabama State University. Specimens were stored at -80°C until time of DNA extraction.

Approximately 10 mg of thawed liver and spleen tissues were placed into Lysing matrix E tubes (MP Biomedicals) containing zirconia and silica beads, 0.5 ml phenol/chloroform/isoamyl alcohol (25:24:1) (TE saturated, pH 8.0) and 0.5 ml of 2× TENS buffer [100 mM Tris–HCl (pH 8.0), 40 mM EDTA, 200 mM NaCl, 2% SDS] (Wan et al., 2011). Tubes were homogenized by mechanical horizontal vortexing in a Mini Beadbeater (BioSpec Products) at speed 40 and time 6, briefly cooled on ice, and centrifuged at 16,000 rpm for 5 min. Supernatants were transferred to 2.0 ml Phase Lock Gel tubes (Invitrogen) containing 0.3 ml of 7.5 M ammonium acetate and equal volumes of chloroform. Tubes were mixed by repeated moderate inverting 10 times and supernatants were transferred into new tubes containing 0.6 volumes of ice cold isopropanol and 3 μl of GlycoBlue Coprecipitant (Life Technologies). After gently inverting several times, samples were incubated at -80°C for 10 min. Following centrifugation at 16,000 rpm for 5 min, isopropanol was decanted and pellets were washed with cold 80% ethanol and allowed to dry for 5 min. Pellets were eluted with 100 μl of TE buffer. DNA was quantified by NanoDrop2000^TM (Thermo Scientific) measuring the absorbance at 260 nm.

V4 and V3-4 hypervariable regions of 16S rRNA gene were amplified for sequencing at RTL Genomics (Research and Testing Laboratory, Lubbock, TX, United States) in two-step, independent reactions using HotStar Taq Master Mix Kit (Qiagen) with universal primers 515F-806R for the V4 region and primer constructs 357wF/785R for the longer, combined V3-4 regions. Primers for the first step were constructed using the fragment-specific forward and reverse primers (515F-806R or 357wF-785R) with the Illumina i5 and i7 sequencing primers added to the 5′-end of each, respectively. Products from the first amplification were added to a second PCR step based on qualitatively determined concentrations (amplicons were run on a 2% ethidium gel, gel bands were scored, and a volume of products was added to the second PCR based on the scores). Primers for the second PCR step were designed using Illumina Nextera PCR primers with 8 bp dual indexes. Each PCR amplification included 9 μl of sterile deionized H₂O, 0.5 μl of 5 μM forward primer, 0.5 μl of 5 mM reverse primer, 1 μl of DNA template, and 14 μl of Taq Master Mix. The negative control was a reaction mixture with no template DNA. PCR reaction conditions included initial denaturation at 95°C for 5 min, then 25 cycles of 94°C for 30 s, annealing at 54°C for 40 s, and extension at 72°C for 1 min, followed by 1 cycle of 72°C for 10 min and 4°C hold. Barcoding PCR reactions were conducted under the same conditions, except with only 10 cycle extensions. Amplification products were visualized with eGels (Life Technologies). Products were then pooled equimolar and each pool was size selected in two rounds using SPRIselect beads (BeckmanCoulter) in a 0.7 ratio for both rounds. Size selected pools were then quantified using Qubit 2.0 fluorometer (Life Technologies) and loaded on an Illumina MiSeq 2x300 flow cell at 10 pM and sequenced.

The sequence data were analyzed using a standard microbial diversity analysis pipeline, which consisted of two major stages, denoising and chimera detection followed by microbial diversity analysis. Denoising was performed using various techniques to remove short sequences, singleton sequences, and noisy reads. Chimera detection was performed using the UCHIME chimera detection software in de novo mode (Edgar et al., 2011). Lastly, remaining sequences were then corrected base-by-base to help remove noise from within each sequence. During the diversity analysis stage, each sample was run through the analysis pipeline to cluster reads into OTUs (at 97% identity) using the UPARSE algorithm (Edgar, 2013), and then globally aligned using the USEARCH global algorithm (Edgar, 2010) against a database of high-quality 16S rRNA gene sequences to determine taxonomic classifications. After OTU selection was performed, a phylogenetic tree was constructed in Newick format from a multiple sequence alignment of OTUs done in MUSCLE (Edgar, 2004a,b) and generated in FastTree (Price et al., 2009, 2010; Shah et al., 2016).

Microbial diversity of cadaver samples was examined from two perspectives using the phyloseq package in R (McMurdie and Holmes, 2013). First, overall richness (i.e., number of distinct nucleic acid sequences present within the microbiome) was expressed as the number of OTUs and was quantified using the Chao1 richness estimator. Secondly, overall microbial diversity, determined by both richness and evenness and the distribution of abundance among distinct taxa, was expressed as Shannon–Wiener species diversity. Measures of microbial diversity were screened for group (region, organ, gender, manner of death, PMI, season, location, weight, and height) differences using an analysis of variance (ANOVA). Multivariate differences among groups were evaluated with permutational multivariate analysis of variance (PERMANOVA) using distance matrices function ADONIS (Oksanen, 2011). For PERMANOVA, ADONIS distances among samples first were calculated using unweighted or weighted UniFrac via the phyloseq package in R (McMurdie and Holmes, 2013), and then an ANOVA-like simulation was conducted to test for group differences. Principal coordinates analysis (PCoA) using unweighted and weighted UniFrac distances and relative abundance bar plots were generated to visualize relationships and differences among groups. All analyses were conducted in R (R Development Core Team, 2010) and all plots were generated using the ggplot2 package (Wickham, 2009).

Bioinformatic characterization of relative abundances and microbial diversity of the thanatomicrobiome was performed through metagenomic analyses in order to determine if there was greater discriminatory power exhibited by 16S rRNA gene V4 versus V3-4 hypervariable region amplicons. Operational taxonomic unit (OTU) data were validated by rarefaction analyses. Rarefaction data confirmed complete coverage until 20,000 sequences to observe all taxa as shown by convergence of vertical asymptotes for all curves (data not shown). Relative abundances of the top most abundant bacteria according to taxonomic hierarchy are shown in Figure 1. The highest percentage of bacteria on the order level was Clostridiales and seven of the top species were Clostridium spp. Furthermore, 95% of samples contained Clostridium spp., whereas six of the seven samples that did not contain Clostridium spp. were from the V3-4 hypervariable region.

Comparison of Chao1 richness estimations was calculated and the V4 hypervariable region had a higher proportion of average calculated estimates than V3-4 (Figure 2). For V4 region amplicons, the average Chao1 richness estimate was 217 species, whereas V3-4 averaged 125 species. Also, ANOVA analysis revealed a statistically significant difference in Chao1 richness between two 16S rRNA gene regions (ANOVA; p < 0.001) (Table 1). In comparisons of gender and manner of death (accident, homicide, natural, suicide, undetermined), statistically significant differences were observed; however, patterns of species richness were statistically independent of time of death. Conversely, no significance was observed in Chao1 richness in the comparison of other variables (e.g., season of death).

In a comparison of the Shannon Diversity within the total microbiome data, both hypervariable regions demonstrated an overall similar profile (Figure 3). The average Shannon Diversity index representing both regions was approximately 2.63 for the V4 region and 2.55 for V3-4; however, significant differences were observed for gender, manner of death, and season of death (ANOVA; p < 0.001, Table 2).

Results of multivariate difference interactions between both 16S rRNA gene regions and other variables resulted in statistical significance in region and location (interaction ADONIS; p = 0.001), which demonstrated that location was the only factor that had a confounding effect on region and that other factors were not confounding results. Results of ADONIS based on weighted UniFrac distances demonstrated significant differences among gender, manner of death, and season (p < 0.001), but not between V4 and V3-4 16S rRNA gene regions (Table 3). Given the fact that significant differences were observed among regions in unweighted but not in weighted UniFrac, the presence or absence of OTUs was more dissimilar than the abundance of OTUs (Table 3).

In order to visualize beta diversity differences between 16S rRNA gene regions, PCoA plots were generated based on unweighted (Figure 4A) and weighted Unifrac distances metrics (Figure 4B). For unweighted UniFrac, there was relatively low variance among two 16S rRNA gene regions, with only 18.10% of variance explained by primary Axis 1 and 7.27% explained by secondary Axis 2. PCoA plots based on unweighted and weighted Unifrac distances were generated faceted by manner of death and season (Figures 5A,B, respectively). Furthermore, samples for the V4 region clustered compactly together more than those for V3-4. For weighted UniFrac, there was more variance among 16S rRNA gene regions compared to unweighted UniFrac PCoA, with 30.07% of the variance explained by primary Axis 1 and 26.92% explained by secondary Axis 2.

The current research defines a new scientific concept, the “Postmortem Clostridium Effect” (PCE), which refers to facultative anaerobic Clostridium spp. that are ubiquitous during human decomposition. There are three dynamics that contribute to Clostridium species’ omnipresence in decaying humans; one factor involves its very fast doubling time. For example, a species found in the present study, C. perfringens, has the most rapid generation time of approximately 7.4 min at optimal temperatures (37–45°C) (Willardsen et al., 1979). The second factor is the bacterium’s proteolytic functions. Clostridium spp. have collagenases that digest native vertebrate collagen fibers which confer the ability to breach colon epithelial surfaces and mucosal layers and transmigrate to proximate tissues (Harrington, 1996; Burcham et al., 2016). The last advantageous putrefactive factor of Clostridium spp. develops via the cessation of the human heart which results in hypoxia (Gevers, 1975; Proskuryakov et al., 2003). A corpse that lacks oxygenated blood can facilitate enteric anaerobic bacteria, naturally found in the colon (e.g., Clostridium spp.), to efficiently and rapidly flourish in the nutrient-rich host. Previous human decomposition studies reported a marked shift from communities dominated by aerobic bacteria to anaerobic at the end of the bloat stage (Hyde et al., 2013). Taken together, these host–microbe factors within decaying biomass lead to the efficient functioning of the PCE.