Cadaver Thanatomicrobiome Signatures: The Ubiquitous Nature of Clostridium Species in Human Decomposition

Human thanatomicrobiome studies have established that an abundant number of putrefactive bacteria within internal organs of decaying bodies are obligate anaerobes, Clostridium spp. These microorganisms have been implicated as etiological agents in potentially life-threatening infections; notwithstanding, the scale and trajectory of these microbes after death have not been elucidated. We performed phylogenetic surveys of thanatomicrobiome signatures of cadavers’ internal organs to compare the microbial diversity between the 16S rRNA gene V4 hypervariable region and V3-4 conjoined regions from livers and spleens of 45 cadavers undergoing forensic microbiological studies. Phylogenetic analyses of 16S rRNA gene sequences revealed that the V4 region had a significantly higher mean Chao1 richness within the total microbiome data. Permutational multivariate analysis of variance statistical tests, based on unweighted UniFrac distances, demonstrated that taxa compositions were significantly different between V4 and V3-4 hypervariable regions (p < 0.001). Of note, we present the first study, using the largest cohort of criminal cases to date, that two hypervariable regions show discriminatory power for human postmortem microbial diversity. In conclusion, here we propose the impact of hypervariable region selection for the 16S rRNA gene in differentiating thanatomicrobiomic profiles to provide empirical data to explain a unique concept, the Postmortem Clostridium Effect.

Prokaryotic 16S rRNA gene amplicon sequences are extensively used in forensic microbiology as reliable biomarkers for taxonomic classification and phylogenetic analysis of the microbiome of death. The thanatomicrobiome, which is defined as microbial succession in decomposing remains (e.g., blood, bone marrow, liver, reproductive organs), can provide evidence concerning interactions between microorganisms and their mammalian hosts. Microbes symbiotically cohabitate with humans during life, but they also participate in the nature and trajectory of decomposition. The host's death introduces chaos in microbial communities as the body becomes an abounding source of nutrients (Mondor et al., 2012). The question then arises, "What hypervariable region(s) of the 16S rRNA gene best profile the shifts that occur in response to the massive proliferation of microbiota after death?" The phenomenon that these signatures are left behind by the corpse provides unique forensic potential to make available trace evidence that can be used in microbial forensics.
Analysis of the very informative 16S rRNA gene is commonly used as a genetic marker for profiling prokaryotic communities (Lane et al., 1985;Woese et al., 1990;Baker et al., 2003;Tringe and Hugenholtz, 2008;Wang and Qian, 2009). Recent postmortem microbiome studies have focused on the 16S rRNA gene Class I, which spans the V4 region (Hyde et al., 2013(Hyde et al., , 2015Damann et al., 2015;Javan et al., 2016b;Metcalf et al., 2016). The V4 hypervariable region is one of the major functional parts of the microbial gene because it encompasses a portion of the "690 hairpin" (Morosyuk et al., 2000;Wimberly et al., 2000) and decoding center (Schluenzen et al., 2000;Morosyuk et al., 2001). The V3 region is categorized in Class II, which is peripheral to the two functional centers of the 16S rRNA gene (Schluenzen et al., 2000;Schuwirth et al., 2005). Studies have shown that V4 is the best region for phylogenetic studies, particularly at the phylum level (Yang et al., 2016).
A key question is, which sub-region (V4 or V3-4) is more effective for phylogenetic studies of the human thanatomicrobiome? To explore the potential to determine cadaver thanatomicrobiome signatures using two hypervariable regions of the 16S rRNA gene, we compared the performance of primers 515F-806R (V4) to 357wF-785R (conjoined V3-4) hypervariable regions (Johnson et al., 2016). We hypothesized that by modulating 16S rRNA gene hypervariable fragment lengths on the Illumina MiSeq platform, the two specified regions would produce dissimilar microbial signatures.
Bioinformatic surveys have shown that hypervariable regions of 16S rRNA gene differ in the detection of sequence diversity; thus, a particular region may function well for ascertaining a spectrum of bacterial taxa whereas a different region may exhibit a distinct degree of taxonomic diversity. The V3-4 amplicon has a read length of twice 250 bp that offers an ideal target for Illumina paired-end sequencing and will provide a suitable framework for V4 region comparisons of the effectiveness of hypervariable region performance. Here, we performed a phylogenetic assessment of species distinctions using V4 versus V3-4 hypervariable regions from postmortem liver and spleen samples from criminal cases. Furthermore, we determined for the first time that fast-growing members of postmortem microbial communities, Clostridium spp., that usually predominate at longer PMIs, also are the most prominent prokaryotes even at shorter time intervals (PMI = 4 h).

Postmortem Sampling of Human Corpses
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.

DNA Extraction of Postmortem Samples
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.

Illumina MiSeq Sequencing
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 Degrees of freedom (df) corresponds to one less than the number of values in the set of means. The p-values are derived from the F-distribution and the significance level is Pr (>F) < 0.001 (shaded cells). using Illumina Nextera PCR primers with 8 bp dual indexes. Each PCR amplification included 9 µl of sterile deionized H 2 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.
FIGURE 3 | Shannon diversity within the total microbiome data. Diversity estimates are colored (red dots, V3-4 regions and blue dots, V4 region) and faceted by 16S rRNA gene hypervariable region group. Mean values (and confidence intervals) in each group are also illustrated. The mean Shannon diversity for both groups is approximately 2.5.

Bioinformatic Analysis
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(Price et al., , 2010Shah 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).

Thanatomicrobiome Sequencing of Postmortem Liver and Spleen Samples
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.

Thanatomicrobiome Alpha and Beta Diversity Analyses
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 Degrees of freedom (df) corresponds to one less than the number of values in the set of means. The p-values are derived from the F-distribution and the significance level is Pr (>F) < 0.001 (shaded cells).
Frontiers in Microbiology | www.frontiersin.org 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.

DISCUSSION
Human antemortem microbiotas are well documented by the Human Microbiome Project (Peterson et al., 2009), for various body locations of healthy individuals; however, currently there is a paucity of knowledge and need for an in-depth interpretation concerning postmortem microbial communities of internal body sites. Our thanatomicrobiome research represents the largest exploratory study examining a cohort of 45 corpses in fresh and bloat stages. Also, the study provides the first extensive catalog of postmortem microbiomes obtained in internal locations analyzed by two different hypervariable regions of the 16S rRNA gene. Approximately 95% of the postmortem liver and spleen profiled in this study involved Clostridium spp. Moreover, the findings revealed that V4 and V3-4 hypervariable testings represent incongruent phylotype diversity and consequently support individual representative assessments of the thanatomicrobiome. For example, studyspecific disparities were observed; Clostridium spp. were not obtained in only one of the V4 region sequences. On the other hand, these species were not obtained in six V3-4 region analyses. Here, we demonstrate that amplicons more sufficient to discriminate Clostridium spp. in postmortem tissue are derived from the V4 hypervariable region. According to previous thanatomicrobiome studies, Clostridium spp. predominated at long PMIs (up to 10 days) (Javan et al., 2016b). However, the current study determined that these Gram-positive, anaerobic extremophiles also predominate at shorter PMIs (4 h).
Our results support Yang et al. (2016) study based on geodesic distances that suggested V4 was the best sub-region for phylogenetic analysis. In the present study, we confirmed that the V4 region, belonging to Class I, had elevated sensitivity for the detection of forensically relevant bacteria; whereas V3 from Class II showed moderate sensitivity. Of particular interest, there was a higher enrichment of the species Methyloferula stellata discovered by targeting V3-4 amplicons compared to V4 (Figure 1). M. stellata is a methanotroph that grows exclusively on methane and methanol (Gill and Landi, 2011;Vorobev et al., 2011). During the bloat stage of human decomposition, methane, along with various odoriferous putrefaction gases, is produced in high abundance by anaerobic fermentation especially emanating from the gastrointestinal tract (Gill and Landi, 2011). Our study is the first to confirm a bacterial taxon that thrives on one of the putrefying gasses produced during decomposition through the use of 16S rRNA gene V3-4 combined regions in human internal body sites. Another very interesting finding was high abundances of three bacteria, Escherichia coli, C. septicum, and a Pseudomonas sp. detected only in female cases using both hypervariable regions (Figure 1).
In the last decade, postmortem microbiology studies have created novel thanatomicrobiome and epinecrotic communities catalogs using expertise in genetics, next-generation sequencing, and bioinformatics. The creation of a Human Postmortem Microbiome Project (HPMP) will facilitate the development of modus operandi used to empower data comparisons obtained from different national and international laboratories. Extension of existing standard operating procedures that cover sampling, processing, sequencing, and analysis will conceive universal standards for microbial analysis to unify the global research community. In addition, the HPMP framework includes research emphases that will provide the scientific community with a hub through which researchers can explore microbial life after death.

Postmortem Clostridium Effect
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.

CONCLUSION
The thanatomicrobiome contributes a substantial function in modulating human decomposition. Studies are needed to elucidate if hypervariable regions are capable to discriminate all bacterial species; therefore, our emphasis was on the phylogenic resolution produced by small subunit surveys to characterize microbial mediators of decay. Conceivably, thanatomicrobiome analysis will be used to build predictive Thanatos models employing the PCE that can further designate the recovery of distinct community types associated with postmortem microbial communities.

ETHICS STATEMENT
The study was approved by the Committee for the Protection of Human Subjects, Alabama State University Institutional Review Board (IRB) number 2016011. Methods were in accordance with relevant guidelines and regulations regarding working with cadavers. Written informed consent was obtained from next-ofkin relatives of the cases.

AVAILABILITY OF DATA AND MATERIAL
All data generated or analyzed during this study are included in this published article and its supplementary information files.

AUTHOR CONTRIBUTIONS
GJ designed the study and collected human corpses. GJ, SF, JM, and TS extracted genomic DNA, PCR, gel electrophoresis the samples. JW performed MiSeq sequencing and data analysis. GJ and SF wrote and edited the article. All authors read and approved the final manuscript.

FUNDING
This work was supported by the National Science Foundation (NSF) grant HRD 1401075.