Characterizing the postmortem human bone microbiome from surface-decomposed remains

Microbial colonization of bone is an important mechanism of postmortem skeletal degradation. However, the types and distributions of bone and tooth colonizing microbes are not well characterized. It is unknown if microbial communities vary in abundance or composition between bone element types, which could help explain differences in human DNA preservation. The goals of the present study were to (1) identify the types of microbes capable of colonizing different human bone types and (2) relate microbial abundances, diversity, and community composition to bone type and human DNA preservation. DNA extracts from 165 bone and tooth samples from three skeletonized individuals were assessed for bacterial loading and microbial community composition and structure. Random forest models were applied to predict operational taxonomic units (OTUs) associated with human DNA concentration. Dominant bacterial bone colonizers were from the phyla Proteobacteria, Actinobacteria, Firmicutes, Bacteroidetes, and Planctomycetes. Eukaryotic bone colonizers were from Ascomycota, Apicomplexa, Annelida, Basidiomycota, and Ciliophora. Bacterial loading was not a significant predictor of human DNA concentration in two out of three individuals. Random forest models were minimally successful in identifying microbes related to human DNA concentration, which were complicated by high variability in community structure between individuals and body regions. This work expands on our understanding of the types of microbes capable of colonizing the postmortem human skeleton and potentially contributing to human skeletal DNA degradation.


Introduction
Human skeletonization occurs following the decomposition of skin and soft tissue, exposing the bones to the surrounding environment [1]. Once the body has progressed to a skeletonized state, teeth and bone become the best materials that can be used for DNA identification. However, while bone is more recalcitrant than soft tissue, it is not stable; it continues to decay over time. With death, bone undergoes decomposition and diagenesis, the postmortem alteration of bone by chemical, physical, and biological factors that result in modification of the original bone material [2]. Time alone is not a good indicator of skeletal DNA preservation [3]. Instead, bone diagenesis and DNA survival are highly dependent on the depositional environment, including microbial activity [4,5], just as the decomposition of all organic resources are influenced by the decomposer community and physicochemical environment [6]. Bone decay mechanistically proceeds via chemical and/or microbial degradation of the organic and inorganic components of bone [7]. Microbes are capable of colonizing and degrading human bone, and microbial DNA is often co-extracted with human DNA, which interferes with downstream processes [8][9][10]. The organic component of bone consists of 90-95% type I collagen (primarily made up of glycine, proline, and hydroxyproline), with minor contributions from other non-collagenous proteins (e.g., osteocalcin, osteoponin, and osteonectin) as well as lipids, mucopolysaccharides, and carbohydrates. The inorganic or mineral component is most similar to hydroxyapatite and consists of calcium, phosphate, carbonate, and to varying degrees sodium [2,11,12,13]. Bone apatite, or bioapatite, can be described as "nature's trash can," as infiltration and substitutions for environmental elements are common [2]. One of the main requirements for lasting preservation via fossilization is a complete shift from a bioapatite composition to a more stable mineral phase, such as fluorinated apatite or fluorine-and carbonate-enriched apatite [2,14].
When not in equilibrium with the surrounding environment, dissolution and recrystallization of bioapatite occurs, allowing microorganisms and enzymes access to the organic phase, resulting in degradation. Similarly, if the organic component degrades by either chemical or biological means, bioapatite becomes more vulnerable to environmental fluctuations and dissolution of the lattice structure is more probable due to new voids in the crystal lattice [2,7,12,[15][16][17]. For example, wet environments exhibit increased rates of DNA degradation because water allows for mineral dissolution and increased hydrolytic damage [18]. The interdependence between the mineral and organic phases of bone supports the idea that greater porosity increases the susceptibility of bone to environmental influences [19,20], both chemical and biological.
Though the reservoir for long-term DNA preservation in bone remains unclear, binding of DNA to bioapatite crystallites seems to be crucial for long-term DNA survival [15]; persistence within osteocytes or other remnant tissues (e.g., from the red bone marrow) may also be possible [21,22]. Gross bone preservation and weathering has been shown to be unrelated to DNA preservation or degradation in some cases [19], while in others, indices of gross preservation are better correlated [20,23,24]. Differences in DNA preservation and degradation by bone type have been observed, though patterns are not consistent between studies (e.g., [9,10,[25][26][27][28]). Whether this has to do with differences in cortical and cancellous bone composition is debated. More porous elements are thought to have increased bacterial presence [20], but increased presence does not necessarily mean increased degradation, as certain microbial taxa may be better adapted to exploiting skeletal material than others.
In archaeology, microbial degradation of bone has been studied primarily through histological research, focusing on regions of microscopic focal destruction [24,[29][30][31][32]. However, culture-based research has shown that collagenase-producing bacteria can use mammalian bone as a substrate (e.g., Alcaligenes pichaudii, Bacillus subtilis, Pseudomonas fluorescens, Clostridium histolyticum) [33]. Others have shown greater macroscopic preservation from archaeological sites with bones lacking culturable collagenase producing bacteria [34]. These observations suggest that bone, and possibly DNA preservation within a bone, may be partially dependent Science Foundation (NSF), to JMD. DOJ and NSF had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Bode Cellmark Forensics provided support in the form of salary for JD, but did not have any additional role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript. The specific roles of these authors are articulated in the 'author contributions' section. The opinions, findings, and conclusions or recommendations expressed in this manuscript are those of the authors and do not necessarily reflect those of the DOJ, NSF or Bode Cellmark Forensics.
on the amount and/or type of microbes colonizing bones. Genera including Pseudomonas, Xanthomonas, Fusarium, and Trichonella have been cultured from bones from diverse archaeological sites [34]. Experimental research has also shown macroscopic destruction phenomena consistent with fungal degraders, specifically the genus Mucor [35], while others have cultured genera from the phylum Ascomycota [34]. Research to date has primarily been limited to culture-based methods, and only a small subset of environmental microbes can be cultivated in the laboratory [36]. Few studies [37,38] have been conducted since the advent of high throughput sequencing technologies, which permit microbial characterization without cultivation. Thus, there is a gap in knowledge regarding the types of microbes capable of colonizing and degrading human bone.
To address this gap, we characterized bacterial and microbial eukaryotic communities from pre-existing skeletal DNA extracts. These extracts were purified for the original purpose of assessing intra-individual variation in human skeletal DNA preservation and because of this, no soil or gut controls were collected for microbial community comparison. The primary goals of the current study were twofold: (1) to identify the types of microbes capable of colonizing different human bone types, and (2) to relate microbial abundances, diversity, and community composition to bone type and human DNA quantity. In doing this we hoped to answer (1) which microbes colonized human bone; (2) what was the intra-individual variation in bone microbial communities (i.e., between elements of a single individual); (3) was there a relationship between microbial communities and human DNA preservation across a limited range of decomposition environments.

Materials and methods
In 2009, three male individuals were placed outside unclothed on the ground surface to decompose naturally at the Anthropology Research Facility (ARF) at the University of Tennessee, Knoxville (UTK) for the original purpose of assessing patterns of human skeletal preservation (S1 Table) [26]. The individuals were donated with consent to the University of Tennessee Forensic Anthropology Center Body Donation Program for the W. M. Bass Donated Skeletal Collection. Because no living human subjects were involved in this research, the Assurance Status of this Project was Exemption Status under section 101(b), paragraph 4, by the University of Tennessee Institutional Review Board. The skeletons were mapped and recovered following complete skeletonization (13 to 23 months), and gently washed with a new toothbrush and tap water at the Forensic Anthropology Center (FAC). Bones were stored at room temperature in paper bags for several months prior to sampling. The same 55 bones and teeth from each individual (total n = 165), which represented all skeletal element types, were selected for sampling (Table 1, S2 Table). All attempts to sample bones from the same site across individuals were made to eliminate sample heterogeneity. Prior to sampling, the external surface of each bone was cleaned by mechanically removing 1 to 2 mm of the outer surface, followed by chemical cleaning via bleach, ethanol, and sterile water. Bones were sampled using a drill and masonry bit at slow speeds; DNA was extracted from sampled bone powder using a complete demineralization protocol [39]. Bone sampling and DNA extraction and analysis were previously described in detail in Mundorff and Davoren [26]. Human DNA quality and quantity were examined to elucidate patterns of DNA preservation by bone type [26]. These remaining skeletal DNA extracts were used in the present study to assess microbial loading via qPCR and microbial community composition and structure using next generation sequencing of the 16S rRNA and 18S rRNA genes.

Microbial and human DNA quantification
As a proxy for total bacterial abundance [40] and colonization of bone, qPCR was used to quantify 16S rRNA gene abundances using the Femto™ Bacterial DNA Quantification Kit (Zymo Research), which targets a 357 bp region of V3-V4. Assays were conducted following manufacturer instructions using a BioRad CFX Connect™ Real-Time PCR Detection System. Samples were quantified in triplicate, while standards were quantified in duplicate, and a minimum of three no-template controls were included in each 96-well plate. Data are presented as gene copy number per gram of bone powder (gene copies gbp -1 ). Human DNA was quantified using the Quantifiler™ kit from Life Technologies (Qf), which targets a 62 bp amplicon; methods and data are reported in [26].

Total DNA quantification
Total DNA was quantified using the Quant-iT™ PicoGreen™ dsDNA Assay Kit (Invitrogen™) using a 200 μL total volume on a 96 well microplate reader. Samples and standards were run in duplicate, with standards ranging from 0 μg mL -1 to 1.0 μg mL -1 . Total DNA concentrations are reported as nanograms per gram of bone powder (ng gbp -1 ).

Percentage cortical content
To calculate the percentage of cortical and cancellous bone in each DNA sample, clinical CT scans of each element were acquired using a Siemens Biograph mCT 64 slice scanner. Scans involved helical acquisition using a 0.6 mm slice thickness, 500 mAs, 120 kV, and bone window with kernel B70s. Data were stored on compact discs and transferred to workstations with image processing software (OsiriX 5.6, Geneva, Switzerland). The DNA sampling site on each element was digitally measured using ImageJ (National Institutes of Health). A macro was created to detect and measure the areas of cortical and cancellous bone (mm) on each CT slice where the sampling site appeared. Measurements of cortical width and height and cancellous width and height were taken separately for each cortical and cancellous bone region for all bone sample sites. Average cortical and cancellous bone width and height measurements were then computed. Due to issues with scan quality, ten of the original 129 samples were removed from the analysis. Percentages of removed cortical and cancellous bone were computed from each DNA sampling site for all elements. Mean percentages of cortical bone composition at each sampling site were divided into seven categories by skeletal element: (1) 80 to 100%, (2) 70 to 79%, (3) 60 to 69%, (4) 50 to 59%, (5) 40 to 49%, (6) 30 to 39%, and (7) 20 to 29%. The first category consists of bones whose sampling sites did not contain any cancellous bone, including the humerus, radius, ulna, femur, and tibia. The second, third, and fourth categories contained only three elements with sampling sites that were composed of over 50% cortical bone. Percentage data were further averaged from each element type across all individuals. The majority of element types revealed consistent measurements between individuals, with standard deviations of 10% or less. Three element types (temporal, occipital, cervical vertebra) exhibited high variability between the three individuals in the relative amount of cortical and cancellous bone removed from the sampling sites.

Next generation sequencing analysis
Total DNA extracts from bone were sent to Hudson Alpha Institute of Biotechnology Genome Services Laboratory (Huntsville, AL) for sequencing of the V3-V4 region of the 16S rRNA gene and V4-V5 of the 18S rRNA gene using 300 PE chemistry on an Illumina MiSeq instrument. Library preparation was performed by Hudson Alpha according to Illumina protocols. Primers included S-D-Bact-0341-b-S-17 and S-D-Bact-0785-a-A-21 for the 16S rRNA gene [41] and 574f and 1132r for the 18S rRNA gene [42]. Raw sequence data is available at NCBI Sequence Read Archive, Accessions PRJNA540930 or SRP194733.
Samples with less than 5,000 reads were removed from analyses, which only affected the 18S rRNA dataset, and remaining samples were rarefied to even depth by the smallest library (16S rRNA minimum library = 48,288 reads; 18S rRNA minimum library = 5,368 reads; rarefied libraries had a mean Good's coverage of 99.2 ± 0.43 and 99.0 ± 0.47, respectively) prior to alpha and beta diversity calculations and ordination (S1 and S2 Figs). Bray-Curtis dissimilarities were computed for all ordinations. Alpha diversity metrics including Inverse Simpson and observed richness were computed using a subsampling approach, in which richness and diversity metrics were computed for a total of 100 iterations, each scaled to even depth.

Sequence quality analysis
Five samples (three reagent blanks and two bone DNA extracts) failed to sequence using 16S rRNA primers, while twenty samples failed to sequence using 18S rRNA primers. Fastqc and multiqc demonstrated high quality reads in the forward direction, with a drop in mean quality Phred scores in the reverse direction at an approximate base pair position of 200 (Phred Score < 25). Following cutadapt and trimmomatic, total 16S rRNA contigs were reduced by 46%. This was further reduced by an additional 14% following further processing in mothur, resulting in a total read loss of 60% (from 37,185,525 to 14,958,201 sequences). This left a total of 14,958,201 sequences, of which 692,709 were unique.
18S rRNA sequences presented an additional challenge; using 300 PE chemistry, forward and reverse reads overlapped by~59 base pairs (bp) (see [48]). Fastqc and multiqc showed a significant reduction in mean base quality in both forward and reverse reads. Forward reads showed a drop in mean quality scores at an approximate position of 250 bp (Phred scores < 25), the same drop in quality was observed in reverse reads at~200 bp. As a consequence, trimming to remove low quality base pairs resulted in a dramatic loss of reads.
Following cutadapt and trimmomatic, total 18S rRNA contigs were reduced by 46%, and after further processing in mothur, sequences were further reduced by 49%, resulting in a total read loss of 95% (from 30,253,173 sequences to 1,518,971). This left a total of 1,518,971 sequences, of which 181,486 were unique. Due to poor read quality, individual A was removed from additional data analysis in phyloseq, resulting in a remaining 7,901 OTUs across 91 samples. Following the removal of samples with less than 5,000 reads, a total of 71 samples remained.

Data analysis
All data analyses, excluding random forests tests, were conducted in R (v.3.4.1). Two-factor analysis of variance tests (ANOVAs) were used to examine differences in log transformed human DNA concentrations by individual and body region (i.e., head, upper torso, arm, hand, lower torso, leg, foot). Assumptions such as normality and homogeneity of variance were tested using D'Agostino's normality test (package = fBasics v. 3042.89) [49] and Levene's test (package = car v. 3.0.2) [50], respectively. Regression analysis was then used to assess the relationship between human DNA concentrations from bone samples and hypothesized predictor variables (i.e., bacterial DNA gene abundances, total DNA concentrations, and percentage cortical content). Human DNA concentrations, bacterial gene abundances, and total DNA concentrations were log-transformed prior to linear regression. Multiple regression analysis was also performed, treating log transformed human DNA as the dependent variable and log transformed bacterial gene abundances, log transformed total DNA, and percentage cortical content as independent variables, including their various interactions. Assumptions including heteroskedasticity, normality, autocorrelation, and multicollinearity were tested using the R package sjstats (v. 0.17.0) [51].
Kruskal-Wallis tests were used to assess statistical significance in alpha diversity metrics, followed by multiple comparisons with false discovery rate (FDR) adjusted p-values. Permutational multivariate analysis of variance tests (PERMANOVAs), applying 999 permutations, were used to assess statistical significance in beta diversity between categorical variables of interest including body region, individual (A, B, and C), human DNA category, and cortical category. These same variables were tested for homogeneity of multivariate dispersion, using 999 permutations. Human DNA category was an arbitrary categorical variable created by dividing a continuous variable, human DNA concentration, by quartiles in each dataset, each quartile defining a category used for factor analysis. Cortical category was established by using the mean percentiles of cortical bone composition at each sampling site as described above [0 (teeth), 1 (80 to 100%), 2 (70 to 79%), 3 (50 to 59%), 4 (40 to 49%), 5 (30 to 39%), 6 (<39%)]. However, because no bones comprised the 60-69% category, this category was eliminated for the purpose of data analysis. In addition, the frontal bone was assigned to the third category rather than the fourth, due to the mean being affected by a single individual. SIMPER, similarity percentages, followed by non-parametric Kruskal-Wallis tests with FDR corrected p-values, were used to determine OTUs significantly contributing to differences between individuals and human DNA category (seq-scripts release v. 1.0) [52]. Random Forest Models were generated using Python (v. 3.5.2) and scikit-learn (v. 0.19.2) [53] to identify OTUs contributing to human DNA preservation. This technique has been frequently applied to forensic microbiome research [54,55] because it has shown success with complex data sets, including microbial abundance data matrices. Moreover, random forest regression has exhibited robustness to overfitting and outputs important features to the model [56], in this case important OTUs. OTUs were merged at the genus level, and all samples were used to generate the model (bacteria, n = 162; microbial eukaryotes, n = 71; combined datasets, n = 71). Data were randomly split into training (3/4) and testing (1/4) sets.

Bacterial and human quantification via qPCR
All no-template controls quantified using the Femto™ Bacterial DNA Quantification Kit were below the limit of detection. The relationship between bacterial gene abundances, used as a proxy for bacterial loading, and human DNA quantities, was not consistent. Despite foot bones having some of the highest human DNA quantities, these also corresponded with high bacterial gene abundances (Fig 1). While bones with high cortical content generally demonstrated lower bacterial infiltration, bacterial gene abundance was not a significant predictor of percent cortical content (adjusted R 2 = -0.03) (Fig 2A). Total DNA was, however, a significant predictor of percent cortical content (p < 0.001, F = 71.43, DF = 1, 33, adjusted R 2 = 0.67); as the percentage of cortical bone decreased, total DNA increased (Fig 2A).
When excluding teeth, human DNA quantities were significantly different by individual (p < 0.001, F = 12.06, DF = 2) and body region (DF = 6, p < 0.01, F = 4.52), with a significant interaction between body region and individual (p < 0.05, F = 1.87, DF = 12). On average, individual B had greater concentrations of human DNA than C or A, with individual A having the lowest concentrations. Therefore, to test the effects of various predictor variables on human DNA recovered from bone, individuals were assessed independently. Bacterial gene When including all predictors (i.e., bacterial gene abundance, total DNA concentration, and percent cortical content) in a single model, the assumption of multicollinearity was not met, indicating that predictor variables were highly correlated.
Eukaryotic communities showed similar distributions in beta diversity compared to bacterial communities. When testing differences between body region, individuals, human DNA quartiles, and cortical content, microbial eukaryotic communities significantly differed by A random forest model was also applied to the microbial eukaryotic dataset to identify OTUs contributing to human DNA preservation. The resulting model was not significant,

Characterizing the postmortem bone microbiome
The postmortem bone microbiome is diverse and variable in the human skeleton two years after death. Excluding Planctomycetes and Saccharibacteria, dominant taxa observed in this study were also shown to dominate human rib samples from twelve individuals that had decomposed at the ARF [37]. Rib samples from the current study most closely resembled dry remains from Damann et al. [37] in phyla-level contributions, but also contained taxa proportions greater than 2% from Verrucomicrobia, Saccharibacteria, Planctomycetes, Chloroflexi, and Chlamydiae (S10 Fig). Discrepancies in observed taxa may be due to differences in sample size and sequencing analysis methodologies. Planctomycetes, a phylum commonly associated with aquatic environments, and Saccharibacteria, a phylum containing multiple environmental taxa have been observed in grave soils [57,58].
Microbial eukaryotic reads succumbed to extensive read removal during sequencing analysis and are less likely to reflect actual eukaryotic community compositions of bone samples. Nevertheless, of the reads that remained, we observed interesting taxa that are worthy of note. Ascomycota, observed in 100% of samples (71 of 71) in the eukaryotic dataset, and Basidiomycota, observed in 55% of samples, were the dominant microbial eukaryotes. This was unsurprising, as these fungal phyla contain multiple saprophytic groups that have previously been observed in association with decomposing carrion [59]. In addition to fungi, multiple phyla of protists were also detected, including Apicomplexa (5 of 71 samples), Ciliophora (49 of 71 samples), and Cercozoa (49 of 71 samples). Protists found in association with bones may be opportunistic, potentially transferred to remains via soil, scavengers, insects, precipitation, and run-off, and may be active fungal and bacterial consumers. For example, the genus Rhogostoma, which was prevalent in samples from individuals A and B, is known to consume both fungal and bacterial species [60]. Similarly, Nematoda were detected, with the majority of sequences belonging to the family Rhabditidae, which contains bacterivorous members, previously observed in decomposition research [55,[61][62][63]. Other bacterivores detected within human bones included Tubulinea, Cercozoa, and Apicomplexa, which have also been found in soils underlying human remains [63]. Cercozoa and other testate amoeba are extremely sensitive to environmental change, and generally decrease in soil with cadaveric inputs [64]. While certain species have responded with positive growth during late stage decomposition (from 1 month to 1 year postmortem) [65], their presence in bones over a year after death likely reflects a shift back to more oligotrophic conditions.
Presence of Deinococcus-Thermus, a phylum well-represented by thermophiles [66], at greater than 1% relative abundance in 6% of samples, is suggestive of a harsh environment.
Bones deposited on the soil surface are exposed to daily and annual temperature contrasts. East Tennessee experiences freezing winter temperatures and temperatures greater than 37˚C in the summer, which can influence moisture availability. As indicated by Reeb et al. [38], bone may provide shelter from harsh environments (i.e., variable temperature, UV). Individual C had greater abundances of Deinococcus-Thermus than B and A (S11 Fig), likely due to the greater duration of exposure to environmental fluctuations, including temperature and precipitation (S9 Fig). The majority of samples with abundances greater than 1% were from the skull including cranial elements and teeth. The cranium is often one of the first anatomical regions to skeletonize during decomposition due to low tissue biomass and high larval presence [67] and likely experiences greater intervals of environmental exposure.

Community differences by individual and anatomical region
Beta diversity analyses showed differences in bone microbial communities, including both prokaryotes and eukaryotes, by individual and body region. This is unsurprising, as there is extensive research on the living human microbiome and the multitude of variables leading to differences in microbial community structure and composition between individuals including life history (e.g., health and diet) [68][69][70]. It is likely that differences in microbial community dispersion influenced PERMANOVA results, and it is possible that individual differences in bone microbial communities would be less evident if a greater number of individuals were included in the analysis. Without other sample types and more individuals to compare, it is difficult to make general conclusions regarding overall bone microbial community similarity. Nevertheless, trends observed across the three individuals examined, suggest that placement duration at the ARF and differences in temperature and precipitation likely contributed to observed differences. In particular, bacterial alpha diversity was lowest in individual A and greatest in C, reflecting differences in exposure duration (S3 and S9 Figs). The impact of soil microbiota is expected to increase overtime with prolonged soil contact [37]. Because of this, we hypothesize that differential rates in skeletonization likely influenced bone microbial composition and structure at any given time point, which may have implications for postmortem interval estimation.
Recently, Pechal et al. [71] showed microbial differentiation by anatomic region (i.e., external sites from the auditory canal, eyes, nose, mouth, umbilicus, and rectum) up to 48 hours after death. Though they speculated that this pattern would likely attenuate with longer postmortem intervals, this has yet to be tested. Here, bone microbial communities retained differences by anatomic location in three individuals with postmortem intervals greater than 1 year. While not tested in this study, micro-environmental differences in soil communities as well as differences in enteric microorganisms and their abilities to compete and persist with soil microorganisms colonizing the body may have contributed to spatial differences observed by anatomic region and between different individuals. Previous research on the human microbiome has shown microbial community uniqueness by individual as well as body site and time [72,73], which has recently gained utility in forensic contexts [74,75].
Nicholson et al. [76] demonstrated that bones in similar environments (i.e., similar pHs and drainage regimes) showed drastic differences in bone preservation. While these parameters could be similar with differences in soil or microbial community composition, there are other factors to consider, beyond the deposition environment, that could potentially impact differential degradation. For example, enteric/putrefactive bacteria have been posited as the primary source of microbial bone degradation in pig remains. Neonatal pig remains demonstrated no evidence of microbial degradation, which researchers hypothesized as being related to the relative sterility of infant guts compared with adults [31]. While the source of bacteria in this study remains unknown, as we have no gut or soil samples prior to placement to track bacterial translocation, we suspect that both soil and gut microbes are able to colonize and aid in bone degradation (e.g., [54]). We have previously demonstrated that human-associated Bacteroides, an obligate anaerobic member of the human gut microbiome, can persist for long time periods in soils exposed to decomposing human remains [57,61], providing evidence that these gut microbes are transferred to the environment and have the potential to colonize bone. The extent to which enteric microorganisms are able to move throughout the body postmortem is likely limited, and distance from the gut may be a crucial factor controlling differences in microbial communities by body region. However, Pechal et al. [71] recently observed an increase in gene abundance associated with bacterial motility during decomposition, so this area of postmortem microbiology merits further study.
Bone microstructure (i.e., the percentage of cortical content) also influenced differences in microbial communities. Communities differed by cortical bone percentage likely due to the presence of greater void space in cancellous bone compared with compact cortical bone, facilitating ease of invasion, especially for incidental taxa or soil contaminants (e.g., potentially Verrucomicrobia). However, this may also be related to nutritive differences; cancellous elements may harbor more labile remnant material such as red marrow [22], while cortical bone may be considered more recalcitrant. This may account for observations in total DNA concentrations and bacterial gene abundances. Bacterial gene abundance was not a significant predictor of human DNA concentration in this study. In cases where bacterial gene abundance did significantly predict human DNA (i.e., individual C), the relationship was positive, indicating that the degree of microbial loading did not negatively impact the pattern of skeletal DNA preservation in remains with environmental exposure up to two years. Rather, the presence of specific taxa likely has a greater impact on skeletal integrity. However, results should be approached with caution, as only three individuals were assessed.
Additionally, presence of both aerobic and anaerobic genera points to the existence of micro-spatial differences within a single bone. This phenomenon is also observed in soils where anaerobic microsites can persist within a well-drained, well-aerated soil. Extracellular polymeric substances were observed surrounding living cells on bison bone at Yellowstone National Park [38]. This highlights the importance of biofilm production in microbial bone colonization. Though microscopy was not performed here to confirm biofilm presence, we hypothesize that biofilm production combined with increased microbial biomass during decomposition plays an important role in the development of micro-spatial differences in oxygen access and respiration strategies.

Microbial taxa associated with skeletal DNA preservation
Random forest models were minimally successful in identifying microbes related to DNA preservation; models were likely complicated by microbial community differences by individual and body region. While bacterial OTUs produced more accurate random forest models than eukaryotic OTUs, the best model resulted when combining both bacterial and eukaryotic data sets, with a Saccharomycetales OTU identified as the most important contributor to the model. Saccharomycetales, commonly associated with the oral microbiome of healthy humans [77], decreased in relative abundance with increased human DNA concentrations in the cranium of individual B. It will be interesting to see the extent to which oral microbes persist in the decomposition environment and impact skeletal degradation as research on the postmortem bone microbiome becomes more prevalent.
Bacterial random forest models were conflated by body region; genera Actinotalea and Paracoccus, showed increased relative abundances with human DNA concentrations in teeth, while Dermacoccaceae demonstrated increased relative abundances in feet. Importantly, increased relative abundances of Clostridium, a genus that contains known collagenase producers [12], were associated with decreased human skeletal DNA concentrations. The foot is the farthest anatomical region from the gut, and interestingly, bones of the feet had some of the highest human DNA concentrations. If the Clostridium present in bones is primarily derived from the gut, then distance from the gut may be an important factor related to human DNA degradation. Of note, increases in relative abundance could just as easily be driven by decreases in other taxa, and so experimental testing with measures of absolute abundance (i.e., qPCR) are needed to validate observed trends [78]. In addition, though predictor taxa could be identified using random forest models, their functional role in DNA degradation, if any, remains unclear. The variation seen by body region and individual may be minimized by increasing the research sample size to include more individuals.

Conclusions
Most of what is known regarding the microbial degradation of bone is from histological research concerning archaeological bone (e.g., [29,32,34,79,80]). The current study used next generation sequencing technologies to provide a survey of bacteria, fungi, and protists potentially capable of bone colonization. Though specific taxa correlated with human DNA preservation using random forest models, the functional role of identified bone microbes remains unknown. Because the target of this study was DNA, which provides information regarding presence rather than activity, it is difficult to discern incidental taxa, i.e., taxa that are present and inactive, from taxa that are actively degrading bone. This is a longstanding challenge in microbial ecology: Linking structure and function. Remnant extracellular DNA of microbial origin is a problem [81], and microbial DNA can bind to hydroxyapatite similar to human DNA [82], further complicating observed differences in community composition and structure.
Additionally, bones were stored at room temperature for several months prior to sampling, likely influencing communities in the bone due to shifting conditions from an outdoor to indoor environment. Though all bones used in this study were stored under the same conditions, we do not know whether bone microbiomes associated with each individual or bone type responded to environmental shifts due to storage conditions in the same way. We have made the assumption that any bias introduced was consistent between samples. Nevertheless, all conclusions should be viewed through the lens of this assumption. The current study presents a first step in characterizing microbial community differences across bone types within and between individuals following skeletonization. No other study, to the authors' knowledge, has characterized the intra-individual variation in the postmortem bone microbiome. Ultimately, this provides a foundation for understanding the postmortem colonization of bone by microbes and its relationship to bone stability and human DNA preservation. Researchers considering the postmortem bone microbiome for PMI estimation, in particular, should be aware of microbial community differences by bone type.  Table. Bacterial alpha diversity metrics including observed richness and diversity (Inverse Simpson). Data was separated by individual, and the mean and standard deviation was computed by body region. Significance levels are represented by asterisks ( � p < 0.05; �� p < 0.01); multiple comparison tests by individual were only conducted using Inverse Simpson indices. Lower case letters in parenthesis refer to body regions. A body region with an exponent corresponding to another body region indicates that those two regions have significantly different diversity indices. (DOCX) S1 Fig. 16S rRNA targeted