Bone biodeterioration—The effect of marine and terrestrial depositional environments on early diagenesis and bone bacterial community

Bacteria play an important role in the degradation of bone material. However, much remains to be learnt about the structure of their communities in degrading bone, and how the depositional environment influences their diversity throughout the exposure period. We genetically profiled the bacterial community in an experimental series of pig bone fragments (femur and humeri) deposited at different well-defined environments in Denmark. The bacterial community in the bone fragments and surrounding depositional environment were studied over one year, and correlated with the bioerosion damage patterns observed microscopically in the bones. We observed that the bacterial communities within the bones were heavily influenced by the local microbial community, and that the general bone microbial diversity increases with time after exposure. We found the presence of several known collagenase producing bacterial groups, and also observed increases in the relative abundance of several of these in bones with tunneling. We anticipate that future analyses using shotgun metagenomics on this and similar datasets will be able to provide insights into mechanisms of microbiome driven bone degradation.


Introduction
Archaeological bones can provide valuable information about past vertebrates, including their diet, evolutionary relationships, demography and even health history (e.g. [1][2][3]). Thus, improving our understanding of the post mortem diagenetic processes that bones undergo and the interaction with the depositional environment is crucial for scientists wishing to better exploit them. Microbial bioerosion is one of the most commonly observed forms of bone degradation [4][5][6]. In this process, the microorganisms attacking the bone may leave characteristic tunnels in the bone microstructure that can be used to visually characterise the type of decay [7]. Two types of tunneling are commonly identified and used for characterisation, namely Wedl and Non-Wedl tunneling. Wedl tunneling is believed to be driven by aquatic microorganisms such as cyanobacteria, or fungal attack; whereas non-Wedl tunneling is believed to be caused by bacteria [7,8]. In theory, microbiological tools should allow those interested in bone bioerosion to obtain more nuanced insights into the diagenesis process. Although prior attempts have proven difficult, due to challenges in cultivating and identifying the communities [9] a few studies successfully managed to demonstrate that a selection of collagenase-producing soil bacteria can grow on bone, including members of the Firmicutes (e.g. Bacillus subtilis, Clostridium sp.) and Proteobacteria (e.g. Pseudomonas sp., Alcaligenes sp.) [10,11]. Fortunately, the past two decades has seen the emergence of DNA sequencing techniques such as metabarcoding and shotgun sequencing, that render it possible to genetically profile entire microbial communities within a sample. Within forensic sciences, several studies have used genetic tools to examine the microbial community associated with decaying cadavers, in order to estimate the post-mortem interval (e.g. [12,13]). However, few studies have focused on the microbial community in articulated bones within cadavers, and fewer still on disarticulated bones. Nevertheless, one study used 16S rRNA gene metabarcoding to characterise the bacterial community in partly skeletonised human ribs [14], and found Proteobacteria to be the most dominant phylum across all samples, followed by Firmicutes, Bacteroidetes, Actinobacteria, Acidobacteria and Chloroflexi. Proteobacteria is a diverse phylum that is associated with any natural environment [15][16][17][18]. Subsequently, Damann and Jans [19] also used metabarcoding to explore the bacterial community of bones degrading within whole cadavers over a period of four years. These authors found that, while Firmicutes was the dominant phylum in bones degrading for less than one year, when bones were allowed to degrade for longer (1-4 years), Proteobacteria dominated the bacterial community. Both studies therefore reported a shift in the bacterial community over time after exposure. Taxonomic characterisation of the microbial composition in bones that have been buried, or left exposed articulated or as fragments, has to the best of our knowledge, not been attempted.
Several metabarcoding studies have also characterised the microbial community in soft tissues of cadavers as well as the surrounding environment [20][21][22][23]. In many of these studies a correlation was shown between the microbial community around a decomposing cadaver, and that within the local sediment, indicating that the sediment microbes could play a role in the degradation of the cadaver [13,24]. For example, Metcalf, Xu et al. [13] examined the microbial composition in the sediment prior to deployment of mouse and human cadavers, and found that approximately 40% of the microbial decomposer community originated from the sediment, although the sediment type itself was not a driving factor in the community development.
In this study, we aimed to contribute to the understanding of microbial degradation of bones during the early stages of diagenesis, and to examine how a marine or coast-near depositional environment influences this community. Unlike humans, bones from animal species typically consumed by humans predominantly enter into the archaeological record defleshed, and thus processed in some way. We therefore used an experimental sample set of defleshed raw, boiled and baked pig bones that were exposed in four different environments associated with submerged prehistoric archaeological sites. More specifically, we aimed to profile the microbial community composition over a time series using 16S metabarcoding as well as look for the presence and possible increases of collagenase producing bacteria within the tunneled bones.

Study environments
The effect of bioerosion on bone exposed in four different depositional environments was studied, with all sites located within the temperate climate zone, in Denmark, Northern Europe (S1 Fig). Two of the depositional environments studied were marine burials placed in different types of sediment (submerged gyttja and submerged sand) 10-30 cm below the sea bed. The other two depositional environments were coastal, of which one was a burial at 30 cm depth in sand a few meters from the sea but permanently above sea water level (terrestrial sand), and the second was exposure on the ground in a tidal zone where the samples were flooded twice a day at high tide (tidal zone). In one environment (the submerged sand) we installed samples during both spring (S) and autumn (A), in order to investigate seasonal variation in the bacterial community.
A thorough environmental analysis of the depositional environments is presented in [25]. A summary of the environmental characteristics of the sites are shown in Table 1. In addition to these, two sediment cores (down to 50 cm) were collected from the two submerged environments for visual description and to measure the total organic content, grain size of the sediment and the sulphate and chloride content of the pore water. The sediment pore water was extracted using a Rhizon SMS (Rhizosphere, NL) attached to vacuum tubes via a syringe set, filtered at 0.45 μm, and the sulphate and chloride content was measured on an ion chromatograph (Thermo Fisher). Sediment samples were dried (at 105˚C) and the total organic content was measured by loss on ignition (after ignition at 500˚C) of the dried samples, where the loss was interpreted as the organic content. The grain size was measured on dried sediment samples (without removal of the organic content) on a Sieve Shaker (Buch & Holm A/S, Denmark) with eleven sieves ranging in mesh size 0.063-2 mm. The shaker was run at maximum speed for 2 min after the sample had been loaded onto the 2 mm sieve. The sediment caught in each sieve was weighted and the accumulated amount of sediment was calculated [26] and visually presented in a logarithmic scale graph.

Experimental setup
Fifty-six bone fragments were used in this study, all originating from either femora or humeri (three whole bones of each type) of domestic swine (Sus scrofa domesticus) obtained from a local butcher. We elected to use pig, for the practical reason that it is a larger mammal from which large quantities of similar bone types can be readily obtained and used without difficulty. As mentioned in the introduction, bones from animal species typically consumed by humans predominantly enter into the archaeological record defleshed, and processed in some way. To mimic this, the bones were either boiled (for 2 h in tap water), baked (in the oven at 200˚C for 3 hours) or left simply defleshed (henceforth termed 'raw'), prior to being cut into fragments of 2 x 2 cm. Details on sample preparation and installation are presented in Eriksen, Table 1. Summary of the environmental characteristics of the four depositional environments. Temperature and pH are reported as the median and range during the exposure period except when only a single measurement was available (indicated as n.m.).

Environment
Terrestrial sand Tidal zone Submerged gyttja Submerged sand Matthiesen et al. [25]. In short, two different systems were used for deployment in the depositional environments. At the marine sites the samples were secured on carbon fiber spears, and at two other sites the samples were sewed into fishing net and either buried or fastened on the ground. When designing the experiment, we deliberately chose to use bone fragments pretreated in this way, as it was consistent with similar material typically found in the Danish archaeological record. However, we acknowledge that in doing so, our samples were disconnected from the gut microbiome that could potentially influence the microbial community [27]. As such we caution against viewing our results as representative of buried material that is not defleshed (e.g. typical local human inhumations). Following installation at the four depositional environments, samples were collected at 4, 14, 28 and 52 weeks (see S2 Fig for details on the sampling periods). At one site (the submerged sand), the study was twice (Autumn and Spring) disrupted after 28 weeks due to unforeseen construction work, thus the time series from this site is not for the full 52 weeks. Fragments were temporarily stored after collection at ca. 5-10˚C in cool boxes (1-4h), until they could be stored frozen (-18˚C). An overview of the samples is presented in Table 2. A few grams of sediment were collected at all the depositional environments in close proximity to the bone samples, at both time of installation, and time of retrieval of the bone samples. We additionally collected sediment for the tidal zone and submerged sand samples at a third period. Unfortunately, the sediment samples from the submerged gyttja environment subsequently defrosted due to a breakdown of the freezer they were stored in, thus we elected not to process them, as we deemed it likely that this would have distorted the true community profile in an unpredictable manner.

Scanning electron microscopy (SEM)
To visualise the potential microstructural damage to the bone fragments caused by microbial and environmental factors, a subset of fourteen of the bone fragments were dried in absolute ethanol, vacuum embedded in epoxy resin and prepared for Scanning electron microscopy (SEM) imaging as described in Turner-Walker [28]. Bones with the longest exposure times, as well as unburied controls, were used for the SEM studies.

DNA sampling and extraction
Post recovery, the bone fragments were visually inspected and photo documented (S1-S3 Figs in Eriksen, Matthiesen et al. [25]). Three control bone fragments (raw, boiled, baked) that had been maintained frozen throughout the experiment following initial pretreatment were also sampled. DNA was extracted from approximately 1000 mg of bone powder sampled (using a rotating drill (Dremel Europe, Bosch Power Tools, NL)) from each bone fragment, as well as 2500 mg of each of the eight sediment control samples ( Table 2). Extraction blanks containing only the reagents from the extraction kit were also incorporated in the experiment. DNA was extracted from the bone powder following Eriksen, Matthiesen et al. [25] Text S3, which in brief involves an EDTA incubation to decalcify the samples prior to extraction using the PowerSoil1 DNA Isolation Kit (Qiagen, US) following the manufacturer's protocol, with minor modifications (400 μL sample from above was added to the PowerBead Tubes). A Tissuelyser (Qiagen, Hilden, Germany) was used for 10 min at 30 Hz instead of a MO BIO Vortex Adaptor). Post extraction, the concentration was measured by Qubit TM (dsDNA high sensitivity DNA assay, Thermo Fisher Scientific) and TapeStation (2200 TapeStation, Agilent Technologies, HS Assay). Table 2. Overview of samples used in the experiment. When no data is available, "nd" is reported. For the Scanning Electron Microscopy (SEM) imaging "None" is reported when the SEM analysis yielded no detectable tunneling, when Wedl attack is observed, "Wed" is reported in the table, while "bac" relates to non-Wedl microbial tunneling.

Bacterial 16S metabarcoding
Prior to bacterial 16S metabarcoding, a subset of the DNA extracts was checked for the presence of PCR inhibitors using a qPCR assay based on a generic mammalian barcoding primer set (16Smam1 and 16Smam2), with PCR conditions described in Taylor [29]. The assay, which looked at how cycle threshold (Ct) values changed across a dilution series of each extract, showed no evidence of any PCR inhibition in the extracts. Subsequently, metabarcoding was principally performed as described in Eriksen, Puetz et al. [30] on the fifty-six bone extracts and the eight sediment extracts, using the 515F and 806R primers that amplify the V4 region of the bacterial ribosomal (rRNA) gene [31,32]. PCR products were visualised on a 2% agarose gel and subsequently pooled at roughly equimolar ratios as determined by gel band strength. Library construction on pooled amplicons was done using the Tagsteady protocol following Carøe and Bohmann [33]. Libraries were sequenced at the Danish National High Throughput Sequencing Center in Copenhagen, Denmark on an Illumina Miseq instrument using 300PE chemistry, with each pool getting ca. 15% of an Illumina MiSeq V3 flowcell (Illumina, US).

Processing sequence data
Unless otherwise specified, default parameters were applied in data processing. Overlapping paired-end Illumina sequencing reads were merged with USEARCH (v. 10.0.240) [34]. Merged reads were demultiplexed with Cutadapt (v1.18) [35] and primers and barcodes were trimmed. After trimming, merged reads had an average length of 252 bp. Quality filtering of trimmed reads was performed with USEARCH using the -fastq_filter command with the -fastq_max-ee_rate option set to 0.01 to remove reads with more than 1 error per 100 bases. Trimmed and filtered amplicon reads were denoised in USEARCH with the unoise3 command, which corrects sequencing errors and removes chimeric sequences. An OTU table with denoised zeroradius OTUs was made by mapping reads to OTUs (hereafter referred to as amplicon sequence variant, ASV). Taxonomy of the ASVs was predicted using SINTAX [36] against the SILVA database (v123) [37]. The resulting ASV table was imported into R [38] for downstream analyses. The ASV table was curated for erroneous ASVs using the LULU algorithm [39] and by removing ASVs that could not be taxonomically classified at Phylum level. Samples with fewer than 1,336 reads were removed. The following packages were applied in R for analyses: LULU, phyloseq [40], vegan [41], DESeq2 [42], PerformanceAnalytics [43], pheatmap [44], and ggplot2 [45]. For alpha diversity estimation, the Chao1 index was calculated on exposed bone samples only, using the "estimate_richness" function from phyloseq. A Pearson correlation analysis was performed using the base R function "cor" with Chao1 indices against sampling points for exposed bones. A Permanova test, using the "adonis" function from vegan package, with a Bray-Curtis distance matrix with pretreatment (boiled, baked, raw), environment, and exposure time as explanatory variables was performed to test what metadata contributed significantly to the variation within the dataset. Ordination with non-metric multidimensional scaling (NMDS from phyloseq function "ordinate") was performed on ASV relative abundance counts, while canonical correspondence analysis (CCA from phyloseq function "ordinate") was performed on a Bray-Curtis distance matrix with Hellinger-transformed ASV counts with the depositional environment as constraining variable.
ASVs that significantly differed in relative abundance between depositional environments were identified using the DESeq function, which identifies log2-fold changes. Of the ASVs varying significantly between depositional environments (p < 0.001), the 50 most abundant ones were extracted for plotting in heatmaps, where both ASVs and samples were hierarchically clustered on the axes.
Lastly, we undertook a preliminary analysis of association between bacterial taxonomy and decomposition of the bone material. Specifically, we looked for changes in abundance of bacteria in the bones deposited at the terrestrial sand, tidal zone and submerged gyttja environments, as these are the three depositional environments in which SEM analysis detected bacterial damage on the bone fragments. Based on the visual analysis of the type of damage observed in the bones from the marine versus the terrestrial environment, we assumed different, environment-specific, types of bacteria may have been active in the biodegradation, and attempted to identify candidate taxa that may be relevant for future study.

Environmental analysis
The results from the organic content, grain size and sulphate:chloride analysis of the two submerged environments are shown in S3 Fig along with a brief visual characterisation. The results show that the upper 18 cm at the submerged gyttja environment consists of sandy gyttja with a high organic content (6-9% total organic matter) whereas deeper in the seabed (18-50 cm) the amount of gyttja decreases and the organic content is only 0-2%. The grain size is fine (<250 μm) to medium (<500 μm). As for the submerged sand environment, the upper 40 cm consists of fine sand with a few plant remains and an organic content below 1%, whereas below 40 cm the organic content increases to up to 3%. The grain size is very fine (<125 μm) and increases to fine (<250 μm) in the deeper layers. The results of the sulphate and chloride analysis suggest that sulphate reduction takes place in both environments. The submerged sand environment shows a linear drop in sulphate content down to 40 cm suggesting sulphate reduction at this depth, with sulphate being supplied from the seawater by diffusion through the top sediment. The submerged gyttja has a reduced sulphate content in the top sediment suggesting sulphate reduction in this part of the sediment profile.

Scanning electron microscopy
Overall, we observed signs of bioerosion (Wedl tunneling) in all the bones from the tidal zone after one year of exposure (S4 Fig), and in the baked bone fragment from the submerged gyttja environment (S5 Fig). In one case, we also observed bioerosion (non-Wedl tunneling) in the boiled bone fragment from the terrestrial sand environment (S7 Fig). We did not detect any attack on the bones from the submerged sand environment (exposed for 28 weeks), in the raw and baked bones from the terrestrial sand or in the raw and boiled bones from the submerged gyttja environment after one year of exposure. We also did not detect bioerosion on the unburied control bones.
Besides bioerosion, we also observed signs of chemical demineralisation and enlarged and ragged osteocyte lacunae, as well as enlarged canaliculi in the baked bone (S5 Fig) from the submerged gyttja environment. However, in general the bones from the two submerged environments were very well preserved, exhibiting only minor demineralisation (e.g. S6 Fig). Some caution is warranted because, as described in the methods section unforeseen construction work on the submerged sand depositional environment prematurely halted the experiment thus bones were only exposed for 28 weeks, as opposed to the full 52 weeks at the other environments. The bones exposed for 52 weeks in the terrestrial sand showed signs of microbial tunneling entering from the periosteal border into a depth of approximately 200 μm (S7 and S8 Figs). The bones were also covered with closely adhering roots (Eriksen, Matthiesen et al. [25], S2-S4 Figs), which were clearly visible in secondary electron images of the unembedded bones, but which could not be observed directly in the cross-sections. However, local demineralisation that likely corresponded to root action, could be seen in the microstructural (SEM) analyses.

Metabarcoding results
From the bacterial metabarcoding data, we obtained a total of 13,399,719 raw unfiltered reads from the fifty-six bone extractions, eight sediment extractions, extraction negative controls and PCR negative control. In total, the negative controls had 0-20 reads, except for one extraction negative control which had 341 reads. The average number of reads per amplicon sample was 209,371 ranging from 2,120 to 1,155,345 (see Table 2). The total number of denoised ASVs was 25,778 prior to LULU curation of erroneous ASVs, of which 7278 were present in the bone fragments. The extraction negative control with the highest number of reads had 30 ASVs assigned, whereas the rest of the negative controls had 0-11 ASVs assigned. Based on the separate clustering in the NMDS ordination (S9 Fig), the negative controls were observed to cluster away from the exposed bone. Thus we believe that the source of contaminant DNA was low level templates in the reagents, rather than inter-extraction contamination, and that these contaminants were masked in the actual samples by the more abundant bacterial DNA present in them [46]. Therefore, they are unlikely to affect our conclusions. In Table 2, an overview of the number of raw reads (unfiltered) and untrimmed ASVs for each sample is presented.

The bacterial community at the different environments.
A Permanova test showed that sample environment (p = 0.001, R 2 = 0.395) contributes most to the total dataset variation, while exposure time (p = 0.001, R 2 = 0.084) also contributes significantly, although to a lesser degree. Pretreatment (boiled and baked) (p = 0.611, R 2 = 0.0200) did not contribute significantly to the variation.
A constrained correspondence analysis (CCA, Fig 1), with environment as a constraining variable shows how the bacterial community within the bones cluster by the depositional environment, with a distinct difference between the marine-associated (submerged sand, submerged gyttja and tidal zone), and the terrestrial environment, and the unburied control bones. The bacterial communities in the bones at the two submerged environments differ from the other environments, in having members of Fusobacteria and Firmicutes as major contributors to the bacterial community. The bone bacterial community at the two submerged environments and the tidal zone both contain members of the Aminicenantes and Spirochaetae. In general, however, regardless of environment, the bone bacterial communities of all samples are dominated by members of Actinobacteria, Bacteroidetes, Proteobacteria, Chloroflexi and Plantomycetes. This can also be seen in the relative abundance plot (S10 Fig). Cyanobacteria are also main contributors to the ordination of the bone bacterial community in all bones (except the three control/unburied bones), including the ones from the terrestrial sand, which is not surprising given they are generally ubiquitous in soil [47] and cyanobacteria are present in all sediment samples (S11 Fig).
The bone bacterial diversity in the control samples (unburied bones) differs from that in the exposed samples, as the community in the latter could clearly be seen to be influenced by the bacterial community at the different depositional environments (Fig 1). The bacterial community in the bones from the terrestrial environment assembles the one in the control sample for the longest period (until week 28) before showing more similarities with the bacterial community found in the environment (Fig 2). The most abundant ASVs present in the unburied bones are only present at lower abundance and decreasing over time in the exposed samples (except for ASV2 belonging to Actinobacteria at the terrestrial sand environment which increases at 28 and 52 weeks) (Fig 3).
The differences in the bone bacterial community composition between depositional environments observed in the Permanova test and CCA plot (Fig 1), were further investigated in the DESeq test of differentially abundant ASVs between environments, and 1011 ASVs were found to be significantly different in their relative abundance between environments (p < 0.001). The top 50 most abundant ASVs that significantly differ in relative abundance between environments are shown in the heatmap (Fig 3). As observed above, the bacterial communities in bones from the submerged environments have many highly abundant ASVs in common. Here is it clear to see that the anaerobic Firmicutes family Lachnospiraceae, and especially the Bacteroidetes family Marinilabiaceae, contribute to the bacterial community in bone fragments at the submerged environments, alongside other anaerobic, sulphate-reducing or marine associated families of the phyla Fusobacteria (Fusobacteriaceae), Proteobacteria (Desulfovibrionaceae) and Thermotogae (Thermotogaceae). For the tidal zone samples, two Proteobacterial families (Helicobacteraceae and Rhodobacteraceae) contribute to the bone bacterial community, as well as several of the same families as seen in the two submerged environments (Desulfovibrionaceae and Marinilabiaceae). Whereas for the terrestrial sand the Proteobacteria families Pseudomonadacea and Xanthomonadaceae, and the Actinobacteria family Nocardiaceae as well as the Verrrucomicrobia family Verrucomicrobiaceae, which contribute to the bone bacterial community.
3.3.2 The bacterial community during the exposure period. Overall, the alpha diversity, as estimated with the Chao1 diversity index calculated on ASVs, of the bacterial community in the bone fragments was observed to generally increase over time (Fig 4; Pearson correlation test). Although this increase with time of the bacterial community diversity is observed in all bones from all environments, the alpha diversity was highest in the bones from the terrestrial sand and tidal zone environments.

The bone degrading bacterial community.
Our SEM imaging analyses revealed that all of the bones exposed for 52 weeks at the tidal zone (ID 8,29,49), and one bone from the submerged gyttja environment (the baked bone exposed for 52 weeks, ID 45) ( Table 2, S4 and S5 Figs) showed similar damage, of a morphology that has previously been associated with Cyanobacteria [8,48]. We therefore analysed this cohort of samples to specifically look for

PLOS ONE
both the presence of, and increases in the relative abundance of Cyanobacteria ASVs. Our analyses indicated that Cyanobacteria ASVs were indeed found in all the bones from the tidal zone, from the submerged gyttja, as well as the sediments in which each was buried. Furthermore, we observed that they had the highest relative abundance in the bones from the tidal zone (Fig 1). We also noted that while Cyanobacteria were absent in the unburied bones, they were present by time of first sampling (4 weeks) in all of the bones from the submerged gyttja and tidal zone.
We next decided to further explore the bone degrading bacterial community within the tunneled bones from the two marine associated environments, by conducting a Pearson correlation coefficient analysis on the relative abundance of ASV's at these two environments. Phyla that significantly correlated with exposure time are shown in S1 Table. We firstly noticed a significant increase in relative abundance occurred with time following exposure, in seventeen phyla from the tidal zone samples, and in twenty-seven phyla from the submerged gyttja environment. Eight of the phyla were shared between the two environments Relative abundance of ASVs, summarised at phylum level, present in the unburied control bone samples and their relative abundance in the sediment and exposed bone fragments (sample IDs on the x-axis are given in Table 2). Plots are subset by both their environment and the exposure time. The three unburied bones exhibit almost no variation and are dominated by ASVs belonging to Proteobacteria, followed by Firmicutes, Bacteroidetes and Actinobacteria. Pretreatment had no significant effect on the variation of the bacterial community (p = 0.611, R 2 = 0.0200). https://doi.org/10.1371/journal.pone.0240512.g002

PLOS ONE
(Nitrospirae, Planctomycetes, Lentisphaerae, Candidate division OP3, Hyd24-12, Chloroflexi, Latescibacteria and Deferribacteres). Notably, none of these phyla are among the most relative abundant in the relevant environmental control samples (S11 Fig). Secondly, we observed that the most abundant ASVs in the bone from the terrestrial sand in which we observed bacterial tunneling (the boiled bone exposed for 52 weeks, ID 33) were the Actinobacteria genera Streptosporangium and Streptomyces (S12 Fig). Several observations suggest that Streptosporangium may be of potential relevance for bone degradation in this environment. Firstly, although it is present at low frequency in both the unburied boiled control bone and the terrestrial sand environment itself (but absent in the other environments), its relative abundance increases in all the bones following exposure in the terrestrial sand environment. Secondly, its relative abundance is highest in the bone sample in which we observed bacterial tunneling.
Thirdly, based on previous studies on collagenase producing bacteria [49,50] we specifically investigated all the bones for the genera Bacillus, Pseudomonas, Bacteroides and Clostridium, along with the Alcaligenaceae family. The Firmicutes genus Bacillus was only present in the sediment at the terrestrial sand environment, as well as at very low abundance in the bones from that environment. However, no correlation between the bone that was tunneled, and the relative abundance of Bacillus was observed. The Proteobacteria Pseudomonas was observed in relatively high abundance in the unburied bones, and in the bones from the terrestrial sand environment. However, it decreased over time in the bones while increased in the sediment. It was only observed, in relatively low abundances, in the bones from the marine associated environments. The relative abundance of both the Bacteroidetes genus Bacteroides and the Firmicutes genus Clostridium within the bones decreased over time during exposure. They were not present in the unburied bones, but observed after 4 weeks of exposure in all bone fragments, but the relative abundance was highest in the bones from the three marine associated environments, only to decrease over time after exposure. The Proteobacteria family Alcaligenaceae were only present in the terrestrial bone and sediment, and were found at highest relative abundance in the one bone fragment that showed tunneling from this environment.  Table 2), nd: no data, none: tunnels absent, Wed = Wedl tunnels.
https://doi.org/10.1371/journal.pone.0240512.g003 Table 1 sums up the environmental conditions at the four depositional environments. While the temperature and pH values are similar, there are fundamental differences in the access of light, oxygen and sulphate at the four sites. All three parameters may influence the activity of the microbial fauna (and thus also the degradation of the bone fragments) by providing energy for their metabolism. The terrestrial sand environment has ample oxygen supply, but no light and only limited amounts of sulphate-rich seawater reaches these bone fragments. The fragments in the tidal zone are exposed to both light, oxygen and seawater. As for the two submerged environments they are not exposed to light or oxygen, but have an ample supply of sulphate and the results indicate a substantial sulphate reduction is ongoing in the most organic sediment layers (S3 Fig). Regarding the difference between the two submerged environments the gyttja environment is a well-known archaeological site [51,52], where the archaeology is preserved in the organic gyttja layers. The site is known to be under erosion [53] and the gyttja layers are relatively thin. The bone fragments were buried where only the upper 18 cm were very organic, below which the gyttja content decreased (S3 Fig). The bones were buried at 5-30 cm depth and thus they may have experienced different organic contents. This intra-site variation is not observed at the submerged sand environment where the upper 40 cm of the sediment are more homogenous. At both submerged environments the uppermost bone fragments will be most prone to degradation, both due to a faster supply of sulphate from the seawater above, and due to the risk of erosion of the sediment.

SEM
We observed bacterial damage on bone fragments exposed for one year, at three out of four depositional environments. In the fourth (the submerged sand) environment, the longest time the bones were exposed was 28 weeks (due to harbour construction work, we lost the rest of the samples at this site), and in those we did not detect any microbial damage (S6 Fig). In the tidal zone fragments and in one fragment (baked) from the submerged gyttja environment, we observed microbial attack (S4 and S5 Figs). Previous studies on bone or tooth submerged in a marine environment have found similar tunneling patterns [54][55][56]. Thus, finding this kind of tunneling in the bones deposited in a marine setting is not in itself surprising. However, the speed in which we observe the damage is interesting. Previous studies have suggested that microbial tunneling does not appear until after four-five years of exposure in a marine environment in the temperate climate zone [57] although in tropical freshwaters tunneling by aquatic microorganisms can appear after only six months [28].
Due to the initial installation design of the experiment, the baked bone fragments at the two submerged environments were buried shallower than the boiled and raw fragments (See Eriksen, Matthiesen et al. [25] Fig 1C). The baked bone fragment from the submerged gyttja environment in which we observe tunneling could have been exposed, or partly exposed, in the water column, potentially making it available for attack by organisms favouring higher oxygen or light levels. This assumption is supported by the aforementioned observations that the whole area is under erosion [53] and thus exposure of the bones during the experimental period is not implausible. This highlights the potential damage cultural heritage artifacts may experience if exposed to the light, oxygenated seawater and the bacteria herein.
At the terrestrial sand environment, we observed irrefutable non-Wedl tunneling in one (boiled) bone fragment after only one year's exposure (S7 Fig). Interestingly, most of the tunneling was confined to a region between 50-200 μm beneath the surface (S8A Fig). This is a common pattern in archaeological bones where the surfaces in contact with the surrounding sediment are relatively un-tunneled compared to the interiors [9]. Perhaps significantly, we observe tunneling that breaches the periosteal surface of the bone at several different locations in this sample (S8B-S8E Fig). The very similar sizes and morphologies of these features suggest a common cause. Although this could be interpreted as simply where advancing bacterial colonies have travelled outwards from the interior and penetrated the surface from below this seems unlikely, considering these tunneling bacteria are generally found to avoid bone tissues close to the bone-soil interface. We therefore suggest that these represent where bacteria have entered the bone from the surrounding sediment and further speculate that the bone degrading bacterial community is to be found among the sediment bacteria.

Bacterial community
Overall, we first explored our dataset by looking for differences in the bacterial community within the bones correlated with environment, exposure time and pretreatment in order to investigate the variation of the bacterial community. We then used the data from the visual analysis as well as previous studies on collagenase producing bacteria to get a better understanding of the bone degrading bacterial community.

Environment.
Our results show that the bone depositional environment contributed significantly to the variation of the bone bacterial community. This indicates that the main contributors to the bacterial community within the bones during degradation were driven by the environment and parameters such as presence or absence of light, oxygen and sulphate, and that the bacterial community within the bones gives an indication of which type of environment the bones have been deposited in. The main contributors to the bone bacterial community at the anoxic environments (submerged gyttja and submerged sand) were obligately anaerobic bacteria (Fusobacteria (Fusobacteriaceae), Firmicutes (Lachnospiraceae) and especially the Bacteroidetes family Marinilabiaceae) (Figs 1 and 3). The two submerged environments were chosen due to their differences in organic content, grain size as well as sulphate: chloride content, to symbolise the degradation within two different kinds of anoxic, sulphate reducing environments. However, the genetic analysis of the bacteria shows that the largest overlap in the bone bacterial communities are between the two marine anoxic environments (submerged sand and submerged gyttja) where light and oxygen are absent. We found that obligate anaerobic bacterial families are main contributors to the bacterial community in the bones at these two environments, and that they share members of Spirochaeta and Aminicenantes with the bones from the tidal zone (Figs 1 and 3).
The bacterial community at the tidal zone is more closely related to the submerged environments than with the terrestrial bacterial community. The marine associated environments all have members of Spirochaeta contributing to the bacterial community. This diverse group of bacteria is widespread in nature, where members occur in a variety of aquatic environments including mud, ponds, lakes, rivers and oceans [58]. Free-living Spirochaeta are facultative anaerobes, which as with Fusobacteria fits with the depositional environments in which the bones were exposed. Aminicenantes are also a very diverse group of bacteria. They have been found across a wide range of salinity and temperature gradients, and it is suggested that they are anaerobic, and capable of fermenting carbohydrates and proteinaceous substrates [59]. Observations from previous studies [14,49] exploring the microbial communities within degrading bones have found that the community predominantly consists of the phyla Proteobacteria, Actinobacteria, Firmicutes and Bacteroidetes, as well to a lesser degree Acidobacteria and Planctomycetes. These studies were done on bones from whole cadavers, however, we find some similarities in the community structure, as across all four environments we find that members of these phyla are the main contributors to the microbial community (Fig 1, S10  Fig). Although we see this similarity in the bacterial community in degrading bones from different environments on phylum level the overall variation of the bacterial community in the bones are environment specific. It seems that environmental parameters such as light, the presence of sulphate as well as oxygen are the driving factors on the composition of the bone bacterial community.
When a bone is exposed to heat, both chemical and physical properties are altered [60]. Prior to obtaining our results we speculated that the heating process would influence the bacterial community within the bones during degradation. This assumption was based on the fact that the degradation of collagen is known to be influenced by heat and water content [61,62] and we incorporated both heating with (boiling) and without (baking) water in our experiment. Heating of the bones is also known to increase the bone porosity, making it easier for the bacteria to access [48]. However, the Permanova test with pretreatment as a constraining variable showed that pretreatment did not contribute significantly to the variation of the bacterial community in the bones. From the SEM analysis, we see bioerosion on bones with both kinds of pretreatments (boiled, baked) as well as the raw bones from the tidal zone after one year of exposure. Our results indicate that the type of environment, rather than the pretreatment, influences presence/absence of bone bioerosion.

Exposure time.
Our results show that the time of exposure contributes significantly to the variation of the bacterial community (p = 0.001, R 2 = 0.088), although less so than the depositional environment (p = 0.001, R 2 = 0.296). Across all environments the diversity of the bacterial community increases over time during the first year of deposition in which we obtained data. After exposure for one year, the largest richness is observed at the tidal zone followed by the terrestrial sand (Fig 4). That the bacterial community in the bones increases during the exposure period is not surprising, as one would presume bacterial activity to increase during decomposition of the bone material.
We speculate that the presence of bacteria in the unburied bones could result from putrefaction happening shortly after death and prior to sampling. As a general tendency, throughout the dataset we observe that the most abundant ASVs present in the unburied bones are only present at very low abundance, if at all, in the exposed bones, and become less abundant during the exposure period (Fig 2). However, the eight most abundant ASVs shared between all three (raw, boiled, baked) unburied bones, are in four cases present in very high abundance in the terrestrial environments, and not at any of the other environments. These four ASVs are all attributed to Pseudomonas, a common Proteobacteria Gram-negative bacterial genus. Members of the Pseudomonas are known to exhibit both lipolytic enzymatic and collagenase abilities and therefore breakdown lipids and collagen [63,64].

Bone degrading community.
Identifying the bone-degraders within the bacterial communities detected in these bone samples is complex, not least when datasets are restricted to 16S metabarcoding (such as ours). We used two approaches to examine the bone degrading bacterial community. Based on the attack patterns observed from the SEM imaging, we firstly looked for known collagenase producing bacteria in the bones with tunneling. We secondly looked at the relative abundance and increases in specific bacterial phyla within the bones from the environments in which we observe damage. This was done based on Metcalf,Xu et al.'s (13) study of the microbial community in degrading cadavers of human and mice, where they searched specifically for decomposers defined as ". . .microbes that differentially increased during decomposition..", and found that 40% of the microbial decomposers were detected in soil prior to the experiment.
Based on the SEM analysis we observed at least two kinds of bioerosion, Wedl and None-Wedl tunneling. The type of damage observed on bones from the marine associated environments (the tidal zone and submerged gyttja) in which we found bone tunneling ( Table 2, S4-S7 Figs) is normally attributed to cyanobacterial attack [48]. With this in mind we screened the samples for the presence and possible increases in the relative abundance of Cyanobacteria. We found cyanobacteria to be highly abundant in the sediment samples from all environments as well as in the bone fragments after 4-14 weeks (not present in the unburied bones), however, the relative abundance of Cyanobacteria was highest in the tidal zone (both sediment and bones). It is also at the tidal zone we observe tunneling in all bones after one year (we did not visually examine the bones for attack earlier than one year, so the attack could have happened before this). The tidal zone environment is characterised by having both the presence of sulphate-rich seawater, oxygen and light, which is optimal conditions for Cyanobacteria [65]. The relative abundance of Cyanobacteria is also high in the one bone from the submerged gyttja environment in which we observed attack, however, not as much as in the tidal zone. Based on this, there seems to be a correlation between the relative abundance of Cyanobacteria and the presence of tunneling in the bones. However, as the relative abundance of Cyanobacteria are already high in the bones from the tidal zone at time of first sampling (after 4 weeks), it could be speculated that this kind of damage could be observed even earlier on the bones.
To further explore the possible bone degrading bacterial community at the two marine associated environments (the tidal zone and submerged gyttja) in which we found bone tunneling ( Table 2, S4-S7 Figs) we conducted a Pearson correlation coefficient analysis on the relative abundance of ASV's that significantly correlated with time. In doing so, we observed that eight of these phyla were shared between the two environments Nitrospirae, Planctomycetes, Lentisphaerae, Candidate_division_OP3, Hyd24-12, Chloroflexi, Latescibacteria and Deferribacteres. Members of all eight phyla have been associated with marine habitat or anoxic metabolic pathways and as such their presence in the bones seems highly associated with the depositional environment and not necessarily biodeterioration [66][67][68]. However, further investigation of these phyla using shotgun genomic data would enable us to explore the possible presence of collagenase activity or other relevant genes associated with bone degradation.
To explore the bone degrading bacteria in the boiled terrestrial sand bone (ID 33), we screened for the most abundant ASVs present after one year, as this is where we observe none-Wedl tunnelling. The genera Streptomyces and Streptosporangium (Actinobacteria) were the two most abundant genera in this bone. Both genera are extremely interesting as they are known to form hyphae-like structures that reassemble fungal hyphae while degrading organic compounds in the soil [69]. In one study, Zaremba-Niedzwiedzka and Andersson [50] found that Actinobacteria accounted for more than 75% of the bacterial reads from a Neanderthal bone, of which the majority were members of the Streptomyces. Interestingly, they did not observe any of the typical DNA damage patterns associated with ancient DNA in the Streptomyces DNA, thus they hypothesised that it must have been derived from living bacteria that were active in the bone degrading process. Both collagenase and protease activity have been identified in members of Streptomyces [50]. We also observed an increase in relative abundance for the Proteobacterial family Alcaligenaceae, although in general at a lower relative abundance than that of Streptomyces and Streptosporangium. However, previous research has shown that members of the Alcaligenaceae are capable of producing collagenase and can use bone as a substrate [10,49]. Although these observations may be interesting in terms of implying which bacteria acts in terms of the bone diagenesis process, given the dynamic nature of many bacterial species' genomes, and the limited taxonomic resolution of metabarcoding data, future studies using shotgun sequenced metagenome data would be needed to clarify the relevance of our observations. Based on previous studies on collagenase producing bacteria we screened all of the bone fragments for the genera Bacillus, Pseudomonas, Bacteroides and Clostridium, along with the Alcaligenaceae family. It is interesting to notice from these results, that there is no relationship between the relative abundance of these genera/family, and tunneling of the bones from the marine associated environments. We acknowledge that screening our dataset for already known collagenase producing bacteria is not the whole story as we speculate that many more groups of bacteria still undescribed could potentially have collagenase active genes. However, looking for already known collagenase producing bacteria is a start which can be elaborated with shotgun data later. Pseudomonas was only observed in relatively low abundance in the bones from the marine associated environments, and although Bacteroides and Clostridium was observed in relatively high abundance within the bones, their presence decreased over time during exposure. We speculate that more than one wave of bone degrading bacteria is present during the early stages of bone diagenesis, with perhaps the first degrading the easily available collagen and lipid sources, and the second degrading the collagen embedded within the mineralised part of the bone structure. This could explain why tunnels are not visible until the second wave is present. Such a change in the microbial community has been explored in various studies attempting to use the microbial community to estimate a time since death interval. Over a longer time period (1-4 years), Damann and Jans [19] also found a change in the bone bacterial community over time. As mentioned earlier, these observations could be interesting in terms of bone degradation, however, while amplicon metabarcoding has its clear advantages being both cost effective, fast and used in a variety of studies making the genetic databases available more exhaustive, the approach still has its disadvantages [32]. With the advancement of full shotgun metagenomic sequencing techniques, the sequence bias experienced with amplicon sequencing and issues of underrepresentation of uncultivable microorganisms in the databases, and therefore difficulties in characterisation, is reduced [70], although shotgun sequencing has recently been found to underrepresent low or high GC phyla with up to two logs [71]. A natural next step with this study would be to look into the functions of the different bacterial strains and examine how the bacterial communities are associated with bone tunneling using a shotgun metagenomics approach.
An understanding of the early bioerosion pathways and the influences from the exposure environment on the composition of the microbial community is of great value in the interpretations of archaeological bone specimens, and in the understanding of the early diagenetic stages of bone bioerosion.

Conclusion
We investigated the bacterial community in bones exposed at four different well documented environments using 16S metabarcoding, and in parallel conducted a visual assessment of the microbial damage pattern using SEM microscopy. We found that the deposition environment contributed significantly to the bacterial community composition within the bones. Not surprisingly, we also observed an increase in the bacterial diversity within the bones over time after exposure, in all four environments, with the community in bones from the tidal zone and terrestrial environments showing the most inter-environmental diversity. Furthermore, we did not observe an effect of pretreatment (raw, boiled, baked) on the bone bacterial community. We observed microbial tunneling in a bone fragment from the terrestrial environment and aquatic microbial attack on all fragments from the tidal zone after one year of exposure which corresponds well with the genetic analysis. Our findings give an insight into the initial bacterial bone diagenetic pathways as well as contribute to the knowledge of the depositional environments' influence on the bacterial community within degrading bones. Caution should be taken when assessing the sediment data, as data was not obtained from all environments at all time points (see Table 2). (TIF) S1 Table. Bacterial phyla in the tidal zone and submerged gyttja. Phyla with a significant correlation with exposure time are given for the two marine associated environments in which we observe tunneling on the bones after one year of exposure. (TIF)