Releasing the microbes from old bones: the effect of different DNA extraction protocols on microbial community profiling

ABSTRACT

DNA-based characterisation of microbial communities can enable those interested in bone diagenesis to address questions relating to the complexity and diversity of said microbial communities. We explored whether biases could be introduced due to differences in the DNA extraction methods used. We investigated the effect of four alternative approaches on the extraction of DNA from bone, in order to examine the resulting effect on the bacterial and fungal OTUs recovered using metabarcoding. We found that the different extraction methods resulted in differences in the microbial OTU profiles generated, both when looking at the fungal and the bacterial communities within the bone. Our results emphasise the need for consistency when working with DNA extraction if comparison of results between different research groups are to be valid. Furthermore, it is clear that future efforts will be needed to determine which methods may provide the most accurate representation of the microbial community in bones.

KEYWORDS:

• Microbial DNA
• bone diagenesis
• metabarcoding
• extraction
• 16S
• ITS

Introduction

There can be many reasons for wanting to study microbes in bones. In fields such as paleo-pathology (Barnes and Thomas 2006; Grif et al. 2012; Kajiya et al. 2010), archaeology (Hollund et al. 2012; Lopez-Costas, Lantes-Suarez, and Cortizas 2016), conservation science (Child 1995), palaeontology (Keenan and Engel 2017; Morales et al. 2017; Morales and Holben 2014) or forensic science (Damann, Williams, and Layton 2015; Handke et al. 2017; Metcalf et al. 2015), knowledge of the microbial composition in bone can be important for answering specific research questions. Furthermore during bone diagenesis (the processes that cause bone alteration during burial (Hedges 2002)), microbes play an important role in the deterioration of the bone structure. Traditionally, bacteria (Hackett 1981; Turner-Walker 2012), fungi (Wedl 1864) or aquatic microorganisms such as cyanobacteria (Bell, Skinner, and Jones 1996; Pesquero, Bell, and Fernandez-Jalvo 2017), have been argued to be responsible for bone destruction.

Characterisation of the microbial diversity associated with mineralised tissue diagenesis has to date mainly been accomplished through histological analysis with light, electron or confocal microscopy (Bell 1990; Garcia-Donas et al. 2017; Le Garff et al. 2017; Morales et al. 2017; Tjelldén 2016; Turner-Walker 2012). The characterisation is based on visual assessment of the tunnelling system. Wedl was the first to describe any form of mineralised tissue bioerosion in 1864, when he noticed alterations on tooth material that had been soaked in well water, and described the changes as being caused by a fungus. These types of tunnels are therefore called Wedl-tunnel after this discovery. Later Hackett (1981) described how bacteria are responsible for microscopic focal destruction (MFD) and identified three types: “linear longitudinal”, “lamellate” and “budded” MFD, based on visual analysis through the microscope. More recently it has been suggested that the tunnels Wedl saw might have been derived from cyanobacteria rather than fungi (Kendall et al. 2018; Turner-Walker 2019). Despite the large amount of research done on microbial bioerosion, there is still a gap in the knowledge about the composition and function of the microbial community in degrading bone (Damann and Jans 2017).

Within the last decade the related fields of metabarcoding and metagenomics have rapidly evolved, and today it is possible to characterise whole microbial communities using DNA based analyses centred around the power of high-throughput sequencing platforms (Caporaso et al. 2011; Siles et al. 2018). DNA metabarcoding in particular, is an efficient and cost-effective approach for profiling the taxonomic composition of multiple microbial species in a single bulk sample (Taberlet et al. 2012). In particular, the 16S rRNA gene is often the phylogenetic marker of choice for profiling bacterial communities (Sinclair et al. 2014; Tringe and Hugenholtz 2008; Ziesemer et al. 2015), while fungal communities are often targeted across the internal transcribed spacer (ITS) regions of the ribosome encoding genes (Begerow et al. 2010; Ihrmark et al. 2012; Schoch et al. 2012). However, while these PCR based methods are incredibly powerful, the quality of the information recovered will always be contingent on one key assumption – that the DNA initially extracted from the samples under study is an optimal representation of the microbial community within the specimen. Whether this is the case has recently been challenged on several grounds. On the one hand, it has been questioned as to whether PCR based methods can provide truly unbiased insights into the microbial community, in comparison to, for example, shotgun metagenomic approaches (Ziesemer et al. 2015). This phenomenon can be driven by both biases inherent in the primer sets themselves, as well as in the context of ancient DNA studies, the differential rates of degradation of different microbial taxa. On the other hand, and of more direct relevance to this study, is the even more basic challenge of recovering unbiased DNA from the samples themselves. This is of particular relevance to samples that might contain both hard to lyse cells – such as those of many fungal and bacterial species – as well as samples such as archaeological bones that are recovered from an archaeological environment, that may contain chemical and mineral components that could inhibit the metabarcoding process (e.g. (Kuske et al. 1998; Lever et al. 2015; Wunderlin et al. 2013; Zhou, Bruns, and Tiedje 1996)).

Although numerous methods have been developed for the extraction of DNA from bone, these are predominantly developed to account for the fact that the majority of DNA within most ancient bones is microbial (Der Sarkissian et al. 2014; Poinar et al. 2006) and not of interest to the goals of most palaeogenomic studies. Thus efforts to improve protocols for extracting endogenous DNA (DNA originating from the bone itself) have been taken to make sure the microbial part is down scaled (Barta, Monroe, and Kemp 2013; Damgaard et al. 2015; Yang et al. 1998). Recently however, some studies have begun to focus on the environmental or pathogenic microbial DNA component in bone, and thus provide suggestions of how the microbial DNA should be extracted to enable community profiling (Margaryan et al. 2018). For example, Der Sarkissian and colleagues (2014) used a shotgun sequencing approach to characterise the microbial community in seven ca. 200 - 13.000 year-old horse bones, after initially digesting their samples in EDTA (Ethylenediaminetetraacetic acid) overnight twice (the pellet obtained after a first EDTA digestion was subjected to an additional overnight EDTA digestion), then applying a silica-based method as described by Rohland and Hofreiter (2007) then modified by Orlando and colleagues (2009) to purify the DNA. EDTA has long been incorporated into protocols for the extraction of DNA from bone, due to its ability to dissolve the hydroxyapatite matrix (Rohland and Hofreiter 2007). For example, Prado et al. (2002) showed that using EDTA as a pre-digestion agent prior to extractions gave up to 50% more endogenous DNA. Given that as with host DNA (Campos et al. 2012) one might expect bacterial cells to be incorporated within the hydroxyapatite matrix, it follows that EDTA pre-digestion may also be critical when targeting nucleic acids from microbes important in the bone degradation process.

As for silica-based nucleic acid purification methods, their appeal lies in that they can improve the purity of DNA recovered from ancient samples by removing inhibitors (e.g. collagen type-I or Maillard products) that may interfere with the DNA extraction and also other downstream processes such as PCR (Kalmar et al. 2000; Rohland and Hofreiter 2007). One further consideration relates to the findings of studies targeting microbial DNA in general, that have argued how incorporation of a bead-beating step during the extraction captures more diversity by mechanically disrupting otherwise difficult to lyse bacterial cells (Santos et al. 2015). However this step can also lead to chimeric products in the sequencing results, due to fragmentation of the DNA molecule (Salonen et al. 2010; Teng et al. 2018; Wintzingerode, Göbel, and Stackebrandt 1997). In light of the above, and a growing academic interest in the microbial communities within degrading bone, we set out to perform an initial exploration on how combinations of the above pre-extraction methods might affect the yield and community composition of microbial DNA in such materials. Given that our intention was to characterise the differences incurred by the different methods, as opposed to the biological relevance of the microbes recovered, and in light of the much more extensive reference databases available to use, we elected to use metabarcoding as opposed to shotgun metagenomics sequencing to generate our insights. Ultimately in doing so, we hoped to provide insights into the potential magnitude of the effect on recovered OTUs caused by different DNA extraction approaches, and thus potentially help facilitate cross-comparisons of metabarcoding studies of ancient bones.

Materials & methods

Sampling

The left fully fused humerus from a female domestic cow (Bos taurus) dated to the eighteenth century was used for the analysis (Collection ID: NHMD-ZMK-226558). It has been stored at room temperature at the Quaternary Collections at the Natural History Museum of Denmark since its excavation from an archaeological site in Copenhagen, Denmark.

Approx. 1000 mg of bone was drilled from the lateral proximal diaphysis close to metaphysis using a Dremel handheld drill. To avoid contamination, the table, surface of the Dremel tool, and the drill head were cleaned with a 5% dilution of commercial bleach (sodium hypochlorite) solution, then rinsed in 96% ethanol before the sample was collected. A new drill bit was used which had been treated with DNA Away^TM (Thermo Fisher Scientific) before use. The bone powder was collected on aluminium foil and carefully transferred to a 1.5 mL Eppendorf® Lobind microcentrifuge tube (VWR). Twelve 50 mg subsamples were made from the powder sample (Figure 1).

Figure 1. Flow diagram illustrating the four different methods tested during extraction of microbial DNA from bone, and the two different amplicon primer sets used. The different methods tested were one without predigestion (termed “raw”), and three treatments incorporating EDTA pre-digestion and filtering through Amicon tubes. These latter 3 methods differed in whether the undigested pellet, and EDTA supernatent were either mixed back together (termed “total”) or kept separately (termed “EDTA” and “pellet”) prior to final purification.

Experimental setup

We explored the effect of four different methods on the performance of extracting nucleic acids from entire microbial communities in bone (Figure 1). For each method, subsequent DNA extractions were done in triplicates to test for consistency. All laboratory work was conducted at the ancient DNA facilities at the Centre for GeoGenetics, University of Copenhagen.

In the first method we extracted DNA from three subsamples (ZM1-3) taken directly from the raw bone powder, which excluded the EDTA predigestion, and are herein referred to as the “raw” method (Figure 1). The inclusion of this approach was to investigate whether the EDTA predigestion was necessary for microbial DNA extraction in bones, or if it was only useful for extracting endogenous nucleic acids.

In the second method, three subsamples (ZM4-6) were initially predigested in 2 mL EDTA at room temperature with rotation for 19 h. The EDTA fraction was subsequently separated from the remaining solids by centrifugation at 13,000 g for three minutes, after which the EDTA supernatant was pipetted out of the tube. This supernatant was subsequently concentrated and purified by centrifugation through an Amicon 4 mL YM-30 filter. After the initial concentration step, the remaining concentrate was further purified using a subsequent wash with 2 mL ddH₂O. This purified supernatant fraction was then returned to the undigested pellet, after which both were extracted together using the PowerSoil® DNA Isolation Kit (MoBio). This method is referred to as the “total” method.

Two additional variations were explored, to better explore the microbial DNA composition of the bone fraction targeted by the EDTA step. Specifically three additional subsamples (ZM7-9) were also initially treated with EDTA as in the second method, however the Amicon purified EDTA fraction was not returned to the pellet. Rather the pellet (ZM7.1, ZM8.1, ZM9.1, referred to as “pellet” method) and purified EDTA fractions (ZM7.2, ZM8.2, ZM9.2, referred to as the “EDTA” method) were processed in separate PowerSoil® DNA Isolations (Figure 1).

DNA extraction

DNA extractions were performed using the PowerSoil® DNA Isolation Kit following the manufacturer's protocol, with minor modifications as follows. After adding solution C1 the samples were incubated at 60°C for 10 min and all centrifuging were done at 13,000 g instead of 10,000 g. Furthermore a Tissue Lyser (Qiagen) was used for 5 min at 30 Hz instead of a MO BIO Vortex Adaptor as originally suggested in the protocol. The PowerSoil® DNA Isolation Kit was chosen principally due to its incorporation of a bead beating step. An extraction blank only containing the solutions from the extraction kit was run parallel with the extraction of the samples. Post extraction, the concentration was measured by Qubit^TM (dsDNA high sensitivity DNA assay, Thermo Fisher Scientific) and Tape Station (2200 TapeStation, Agilent Technologies, HS Assay) but showed DNA concentrations under the detection limit (<2 ng). The bacterial communities were then PCR metabarcoded using primers that target the V4 region on the bacterial 16S rRNA gene, and the fungal ITS1 region (details next section).

Metabarcoding

Metabarcoding was performed on each DNA extract in triplicate reactions, and with PCR triplicates resulting in a total of 36 amplicons (3 × 3x4) for each primer set used (Figure 1). Extraction and PCR negatives were also run with both primer sets using the same conditions as for the samples.

Bacterial metabarcoding was performed using the 515F and 806R primers to amplify the V4 region of the bacterial 16S ribosomal RNA (rRNA) gene (S1 Table) (Gilbert, Jansson, and Knight 2014; Ziesemer et al. 2015). The primers were 5’ nucleotide tagged (Binladen et al. 2007). Tags were 6 nucleotides long, differed by at least 3 nucleotides, and had 1 or 2 nucleotides added to their 5’ ends to increase complexity on the Illumina flowcell. For fungal metabarcording the universal genetic barcode, ITS1 (199 / 260) was used to amplify the ITS region (Ihrmark et al. 2012). These primers amplify a ca. 350 bp long fragment and the tags used are shown in S2 Table. Prior to metabarcoding, a subset of the DNA extracts were screened using qPCR to check for PCR inhibition. Subsequent metabarcoding PCRs were performed in 25 µL reactions, with each reaction containing 2 µL of template DNA, 2.5 U AmpliTaq Gold polymerase (Applied Biosystems, Foster City, CA), 1× AmpliTaq Gold buffer, 2.5 mM MgCl₂, 0.25 mM of each dNTP, 0.4 mg/mL bovine serum albumin, 0.4 µM tagged forward primer and 0.4 µM tagged reverse primer. Thermocycling conditions were as follows (16SV4): initial denaturation at 94^◦C for 3 min, followed by 33 cycles of denaturation at 94^◦C for 45s, annealing at 50^◦C for 60s and elongation at 72^◦C for 30s, followed by a final elongation at 72^◦C for 10 min. The PCR conditions for ITS were as follows: initial denaturation at 95^◦C for 5 min, followed by 38 cycles of denaturation at 95^◦C for 30s, annealing at 55^◦C for 30s and elongation at 72^◦C for 60s, followed by a final elongation at 72^◦C for 5 min.

Illumina library preparation

After visualisation on a 2% agarose gel, amplicons were mixed at equimolar concentrations in either pools of 13 (for 16SV4) or 19–23 (for ITS1) samples. Concentration was measured using a 2200 TapeStation. PCR negative controls were also added to each of the pooled libraries. As the gel electrophoresis showed no sign of DNA in the extraction negatives, no further work was conducted on them. Amplicon pools were purified using Agencourt AMPure XP (Beckman Coulter, Inc, Brea, California) bead purification, following the manufacturer’s protocol, prior to being converted into Illumina compatible libraries using NEBNext library kit E6070L (New England Biolabs, Ipswich, Massachusetts) and the blunt-ended library adapters originally described by Meyer and Kircher (2010). Each library was PCR indexed with Illumina single index primers, after an initial qPCR to determine the appropriate optimal number of PCR cycles required. The libraries were amplified in 100 μL reactions, with each reaction containing 20 μL of template DNA, 10 U AmpliTaq Gold polymerase (Applied Biosystems, Foster City, CA), 1× AmpliTaq Gold buffer, 2.5 mM MgCl2, 0.2 mM of each dNTP, 0.2 µM IS4 forward primer and 0.2 µM reverse primer with sample specific 6 bp index. The PCR conditions were as follows: initial denaturation at 95°C for 12 min, followed by 12–16 cycles of denaturation at 95°C for 20s, annealing at 60°C for 30s, and elongation at 72°C for 40s, followed by a final elongation at 72°C for 5 min. Following amplification, libraries were purified using Agencourt AMPure XP bead purification, following the manufacturer’s protocol, and eluted in 50 µL of EB buffer (Qiagen, Hilden, Germany). Final concentrations were measured on a Qubit^TM and on the 2200 TapeStation, and then pooled for sequencing on an Illumina MiSeq platform using 250PE chemistry at the National High-throughput DNA Sequencing Centre at the Natural History Museum of Denmark. Each pool was sequenced over ca. 15% of a MiSeq lane.

Processing sequencing data

Raw paired-end (PE) Illumina MiSeq reads for the 16S rRNA and ITS1 amplicons were treated using largely the same approach with USEARCH (v10.0.240) (Edgar 2010). Following merging of PE reads with USEARCH, trimming of sequencing barcodes, adapters and primers was performed with Cutadapt (v1.8.3) (Martin 2011) with the requirement that sequences to be trimmed were found in reads with more than two thirds of their length. Merged and trimmed reads were quality filtered with maximum number of expected errors in each read of 0.5. Quality filtered reads were dereplicated and subsequently clustered to Operational Taxonomic Units (OTU) in USEARCH with the -cluster_otus command with default parameters (minimum 97% similarity between OTUs), where chimeras were also filtered. Taxonomic classification of OTUs was performed with SINTAX (Edgar 2016) against the SILVA (v123) and UNITE (v7.1 “utax” version) databases for 16S and ITS1 regions, respectively, with -sintax_cutoff option set to 0.8. ITS1 OTUs that were not classified with SINTAX were manually curated for taxonomy using BLASTN (Altschul et al. 1990) against the UNITE database (minimum 90% identity threshold). PCR triplicates with less than 400 reads were removed from both datasets prior to statistical analyses. This was the case for four samples in the bacterial dataset which lost one PCR triplicate each during filtering, and for two samples in the fungal dataset which lost two PCR triplicates. One sample lost one PCR triplicate as it contained 95% plant taxa, and was therefore removed from the fungal OTU table. Unclassified and Plantae OTUs were removed prior to construction of the Venn diagrams, which were created using MetaCoMET based on presence or absence of OTUs in the different methods (Wang et al. 2016).

Microbial taxonomic composition and statistical analysis

Statistical analyses were performed using R version 3.4.3 (R Core Team 2014). Relative abundance plots of the top 20 most abundant taxonomic groups at the phylum and genus level were generated for each method using the RAM package. Alpha diversity measures (Observed, Chao1, Shannon and Simpson) of bacterial and fungal communities were calculated with the R package “phyloseq” (McMurdie and Holmes 2013). These alpha diversity estimates were calculated on OTU data that was not subsampled to equal sampling depth, since a potential loss in sensitivity, as discussed elsewhere (McMurdie and Holmes 2013) could not be tolerated. Differences in richness and evenness were assessed using Kruskal–Wallis followed by pairwise tests among methods using a Wilcoxon rank-sum test followed by a Benjamini-Hochberg (FDR) correction for multiple testing (Benjamini and Hochberg 1995). Pairwise comparisons of microbial communities between each of the methods were assessed with the Metacoder package in R (p-value <0.05) (Foster, Sharpton, and Grünwald 2017). The heat trees generated from Metacoder allow for the visualisation of differences between the median proportion of reads for samples from different methods using a Wilcoxon rank-sum test followed by a Benjamini-Hochberg (FDR) correction for multiple testing (Benjamini and Hochberg 1995). In order to enable between-sample beta diversity statistics with Bray–Curtis dissimilarity, samples were subsampled randomly to 5000 reads with the USEARCH -otutab_norm function. The degree of dissimilarities in the bacterial and fungal community composition were quantified among the different methods by a principal coordinate analysis (PCoA) with the Bray Curtis dissimilarity matrix using the phyloseq package in R. A permutational multivariate analysis of variance on Bray–Curtis dissimilarities (Anderson 2001) (adonis; 4999 permutations), as implemented in the R-package “vegan”, was used to test if there were any significant effects of sample pretreatment on bone bacterial and fungal community structure.

Results

Sequence quality

A total of 496,896 reads were obtained from the 35 bacterial metabarcoding reactions sequenced (and three PCR negative controls). The average number of reads per amplicon was 14,197, ranging between 404 and 76,752 number of reads per sample (S4 Table). The total number of clustered OTUs was 375 (Table 1). We also obtained 707,539 reads for the 33 fungal metabarcoding reactions sequences (Table 1). The average number of reads per sample was 21,441, ranging from 2379 to 96,799. The total number of clustered OTUs was 67 at 97% identity threshold, however the majority of the OTUs (42) were initially unclassified, and had to be manually assigned (by BLASTN the OTUs against the online UNITE database) to taxonomic classification (14 of the 42 OTUs was successfully classified at 90% identity threshold, S3 Table). Both the PCR negative control and extraction blank showed no products on the TapeStation, and DNA concentrations were below the Qubit detection limits. The PCR negative controls yielded no reads.

Table 1. Sequencing results. The bacterial / fungal number of merged and quality filtered reads, number of samples and number OTUs per different method and as total for all methods are listed in the table. Some OTUs are present in more than one method. The number of OTUs stated for “All method” are the total no. of unique OTUs.

Method	No. of reads	No. of samples	No. of OTUs
“raw”	5,994 / 197,617	9 / 9	271 / 35
“total”	86,829 / 204,213	8 / 9	147 / 21
“pellet”	127,277 / 245,186	9 / 9	144 / 25
“EDTA”	276,796 / 60,423	9 / 6	107 / 11
All methods	496,896 / 707,539	35 / 33	375 / 67

The different treatments tested were as follows: “raw” bone powder, which excluded the EDTA pre-digestion, “total” which contained bone digested with EDTA, “pellet” was treated with EDTA but the EDTA fraction was removed and “EDTA” which contained the EDTA leftover from the “pellet.”

Bone bacterial communities

Taxonomical composition of the bacterial community at the phylum and genus levels was consistent across the replicate extractions, and triplicate PCRs per extract, for the “total” (ZM4.1-3; ZM5.1-3; ZM6.2-3) and “raw” methods (ZM1.2-3; ZM2.2-3; ZM3.2-3) (S1a and S1b Figures). The rarefaction curves (S3a and S3b Figures) showed that the number of reads obtained from each of the samples covered most of the bacterial community diversity within each of the methods, despite the lower number of reads obtained for “raw”. Overall, alpha diversity was significantly higher in “raw” compared to all of the three methods using EDTA digestion for both richness (observed and Chao1 all FDR < 0.05; Figure 2a) and evenness (Shannon and Simpson all FDR < 0.05; Figure 2a). “EDTA” showed the largest within-method variation in both Shannon and Simpson evenness measures.

Figure 2. Alpha diversity of the a) bacterial 16SV4 region and b) fungal ITS region of the four methods (non-normalized and minimum cut-off of 400 reads).

The “raw” method showed the greatest number of bacterial differences at all taxonomic levels when compared to each of the other three methods that incorporated an EDTA predigestion (top row in Figure 3a). In comparison, all three methods using EDTA predigestion had few differences among bacterial community composition (Figure 3a).

Figure 3. Heat trees showing the pairwise comparisons of a) bacterial and b) fungal communities between the different methods of extracting microbes in bone. The gray tree on the lower left functions as a key for the smaller unlabelled trees. The colour of each taxon represents the log-2 ratio of median proportions of reads observed for each extraction method. Only differentially abundant taxa are coloured, determined using a Wilcoxon rank-sum test followed by a Benjamini-Hochberg (FDR) correction for multiple comparisons. Taxa coloured green are enriched in the method shown in the row and those coloured brown are enriched in the methods shown in the column.

At the phylum level, the bacterial community was dominated by Firmicutes in all methods using EDTA predigestion, with the exception of one sample in the “pellet” method, and four samples in the “EDTA” method, where Proteobacteria and Bacteroidetes dominated (Figure 3a and S1a Fig.). Pairwise comparisons between each of the three EDTA predigestion methods showed relatively few differences at the phylum level; Proteobacteria were more abundant in method “EDTA” whereas, “pellet” and “total” methods showed more Actinobacteria (Figure 3a). The taxonomic composition of the extractions that were performed directly on the bone powder (“raw”), were strikingly different from those that included an EDTA predigestion, in which Proteobacteria, Actinobacteria and Bacteroidetes composed a large proportion of the phyla (Figure 3a and S1a Fig). Pairwise method comparisons confirmed that the “raw” method was more abundant in the above three phyla compare to the EDTA predigestion treatments, which were themselves consistently and significantly more abundant in Firmicutes (Figure 3a).

At the genus level, the most abundant taxa were Paenibacillus and Psychrobacillus within method “pellet”, whereas Bacillus and Aneurinibacillus were the most abundant among samples within the “total” method (S1b Fig). Pairwise comparisons between each of the three methods using EDTA predigestion showed that “pellet” and “total” had higher overall proportions of Paenibacillus and Psychrobacillus than method “EDTA” (Figure 3a). Additionally, Micromonospora was higher within method “pellet” than all other methods, whereas Cohnella and Bacillus were more abundant in method “total” (Figure 3a). The “raw” method samples had over 50% of unclassified taxa below the family level, otherwise Streptosporangium was the most dominant genus.

The Venn diagram (Figure 4a) supports the statistical data in showing that “raw” contains the highest number of unique OTUs (141 OTUs, 37.6%), whereas “total” has 26 unique OTUs (6.9%), “pellet” has 21 (5.6%) and “EDTA” has 17 (4.5%). It is worth noting however that only 35 of the 375 OTUs were shared between the methods. Furthermore, only 17 OTUs (4.5%) were shared between the three methods using EDTA digestion (“pellet”, “EDTA”, “total”).

The different methods had a significant effect on the bacterial community structure (adonis; p = 0.001). The Bray–Curtis PCoA plot showed three distinct clusters, which were confirmed with a gap statistical analyses (Figure 5a). The majority of the samples from methods “EDTA”, “pellet” and “total” clustered together, with the exception of the four “EDTA” and one “pellet” sample that were identified to have differential taxonomic abundance above. All of the samples from the “raw” method clustered tightly together without any overlap with samples from other methods.

Bone fungal community

The rarefaction curve showed that the number of reads obtained from each of the samples covered most of the fungal community diversity within each of the methods. Overall, alpha diversity was significantly higher in the “raw” compared to the other three methods for richness (observed and Chao1 all FDR < 0.05; Figure 2b) and evenness (Shannon and Simpson all FDR < 0.05; Figure 2b), and the “EDTA” method had significantly lower richness (observed and Chao1 all FDR < 0.05; Figure 2b) and evenness (Shannon and Simpson all FDR < 0.05; Figure 2b) in comparison to all other methods. The “pellet” method showed the largest within method variation in all alpha diversity measures (Figure 2b).

At the phylum level, the fungal community was largely dominated by Ascomycota in all methods (S2a Fig), with the exception of two of the “EDTA” samples, which were dominated by Basidiomycota. Pairwise method comparisons confirmed that “EDTA” samples had more Basidiomycota and fewer Ascomycota when compared to the other three methods (Figure 3b). Taxonomically, there were no obvious differences identified among “raw”, “total” and “pellet” methods at the phylum level (Figure 3b).

At genus level the most abundant genera were Pseudeurotium and Aspergillus in all of the methods except “EDTA” (S2b Figure). Within the “EDTA” method, each sample was dominated by a different genus which, collectively, were not present in any of the other methods. Similar to the bacterial dataset, the “raw” method showed the highest abundance of unclassified fungal groups.

The Venn diagram (Figure 4b) supports the statistical data showing that “EDTA” is the method that shares the fewest OTUs with any of the other methods. Again it is worth noting that the number of shared OTUs among methods are one OTU out of the total 76 OTUs. This means that less than 2% of the total amount of OTUs are shared between the methods. No OTUs were shared between the three methods using EDTA predigestion (“pellet”, “EDTA”, “total”).

Figure 4. Venn diagram showing the shared and method specific OTUs for the a) bacterial and b) fungal communities.

The different methods had a significant effect on the fungal community structure (adonis; p = 0.001). The Bray–Curtis PCoA plot showed that “EDTA” formed a tight cluster that was separate from the other three methods, and was responsible for most of the variation along the first principal coordinate axis (Figure 5b). Samples from the “pellet”, “total” and “raw” methods had overlapping clusters (Figure 5b), corroborating the results from the differential abundance analysis visualised in the heat trees (Figure 3b).

Figure 5. PCoA plot based on Bray-Curtis dissimilarities of the rarefied and log+1 transformed a) bacterial 16SV4 region and b) fungal ITS region.

Discussion

The principal aim of this study was to compare the effect of four different approaches for extracting DNA from archaeological bone remains, on the resulting microbial metabarcoding community profile recovered. Specifically, we aimed to examine the degree to which the relative abundance and diversity of the microbes changed according to each method, and not what the biological relevance was of the microbes recovered. In this regard, we draw the reader’s attention to the findings of several key previous studies in the efficacy of methods for profiling microbial communities from both modern and ancient materials (Eisenhofer et al. 2019; Salter et al. 2014; Ziesemer et al. 2015). And in particular that should biologically meaningful information be required, then as long as suitable reference databases are available to compare against, shotgun metagenomics approaches are likely less inherently biased, thus more likely to provide accurate community reconstructions. PCR based amplification methods are also in general limited by the design of the primers which in this case are not able to amplify fragments shorter than the amplicon length (here ca 290 bp for bacteria and ca 350 bp for fungal), and thus any very degraded or fragmented DNA will not have been PCR amplified and sequenced. We also note that the bead-beating step used in all extractions, while shown in previous studies to give a more accurate representation of the microbial community (Santos et al. 2017; Teng et al. 2018), could fragment the DNA in the samples. However as it was used in all methods tested, its effect should be equal in all cases.

Nevertheless, although our data was based on metabarcoding, it does not change the fact that major differences were observed, and thus support the key point that there is a large effect of methods on the OTUs recovered. Specifically, overall our findings showed that there was a large difference in the microbial community composition when using different approaches prior to microbial DNA extraction of bone material. In the bacterial dataset we observed that when samples were not subjected to an initial EDTA digestion (“raw” method) - and to a lesser extent when only the supernatant was extracted after EDTA digestion (“EDTA” method) - the largest difference from the other methods was shown. In the fungal dataset we saw the largest difference from the other methods when only the supernatant was extracted after EDTA digestion (“EDTA” method).

The bacterial community

With regards to the bacterial community data, it was notable that when the bone was not initially predigested with EDTA to dissolve its inorganic component (“raw”), its community structure was most diverse. However given that the majority (>50%) of this diverse community could not successfully be classified, and was significantly different from the taxa found using the other methods, it seems that most of this community comprises of taxa that are rare. Although there is a higher overall diversity and consistency in the PCR triplicates using the raw method, the dominating bacterial phyla (Figure S1a) found have previously been associated with environmental contaminants such as soil (Margesin et al. 2017; Zaremba-Niedzwiedzka and Andersson 2013). Furthermore, given that quantitative real-time PCR analysis of the extracts indicated that the “raw” extracts contained up to 1000-fold less bacterial DNA than in the extracts from the three EDTA based methods (data not shown), we speculate that their signal is subsequently lost among the flood of more common microbial taxa recovered when fully digesting the bone with EDTA, hence why their overall diversity appears lower.

Our findings show very little overlap of OTUs between the different methods. However, when observing the taxonomic variation among the three methods using EDTA digestion, we do not see any significant differences, and they cluster together on the PCoA plot. The difference between the pre-digestion methods seen in both the microbial community (4.5% shared OTUs) and the fungal community (no shared OTUs) could also be explained by the findings of an earlier study by Campos et al. (2012), which argues that the pellet formed after treatment with EDTA mainly consists of the organic part of the bone (collagen), whereas the EDTA soluble supernatant contains the hydroxyapatite. Given that the principal difference between the three EDTA related methods was whether the EDTA and pellet fractions were kept separate or recombined, the observed differences may simply relate to the microbes in the inorganic vs organic components of the bone. Again we speculate whether filtering the supernatant through the Amicon filters may wash away some DNA that is therefore not added back to the “total” method, nor retained in the “EDTA” method. This could explain the lack of OTU overlap found between methods in both the bacterial and fungal dataset.

However, in this regard, a natural intuition would be that as the “total” method tested combines both the EDTA and pellet fraction, this should encompass the variation seen within the two separate pellet and EDTA fractions, which appears not the case. In our study we speculate again that perhaps the observed diversity differences are driven by rare microbial taxa, that simply vary in location throughout the bone, whereas those most common drive the overall clustering of the three EDTA related methods. It could also be that the very low amount of DNA obtained in all of the extractions from all four methods is limiting the detection of diversity and not the diversity itself.

The four most common bacterial phyla found were Firmicutes, Proteobacteria, Actinobacteria and Bacteroidetes. McIntyre et al. (2017) observed in their study that there seemed to be a sequencing bias towards Firmicutes, Proteobacteria and Actinobacteria most likely because they are better represented in databases than others. However other studies have found that Firmicutes and Bacteroidetes are the dominating phyla related to degradation of bones together with Proteobacteria and Actinobacteria (Damann and Jans 2017). Metcalf et al. (2013) also found Firmicutes dominated the initial bacterial community associated with degrading mouse cadavers, however the study focussed on soft tissue decomposition and not bone material. Actinobacteria are also known to be present in sediments where they play a crucial role in the degradation of biomaterials (Goodfellow and Williams 1983; Margesin et al. 2017). That we find these phyla in our bone sample is therefore not surprising. However, as we are unaware of any genetic studies on the bacterial communities that naturally inhabit bones stored in museum facilities, we cannot rule out that they derive from this environment. Although we feel this is unlikely given that the water content of museum bones is low, which is hardly conducive to bacterial growth. It could also be argued that some of the main bacterial genus we observe in our treatments are derived from contamination as e.g. Paenibacillus is observed in three out of four methods are a known laboratory contaminant (Eisenhofer et al. 2019). As we only included an extraction blank after the EDTA pre-digestion step we cannot completely rule out that this is the case. However, in this regard, we wish to highlight that if it indeed is a lab contaminant, then it would be expected to show up in all treatments as the extraction kit used was from the same batch and the extractions was carried out on the same day. Also we did not observe any contamination in our extraction blank both right after extraction (Qubit) and after PCR amplification (gel visualisation).

Whether these four different approaches give the correct picture of what is actually present in the bone of bacterial and fungal DNA could therefore be debated. However, as a predigestion step is often used in studies that aim to examine the microbial community in bone using EDTA (Barnes and Thomas 2006; Der Sarkissian et al. 2014), exploring the effect of these pretreatment methods is clearly relevant.

The fungal community

In general, the fungal dataset had a very low level of OTUs (67) compared to the bacterial dataset (375). This could simply be due to a less diverse bone fungal community, however the lack of taxonomic diversity in the comparative databases could also lead to an underestimation of the fungal taxonomical identification. The fungal dataset show very few overlapping OTUs (<2%) (Figure 4b) between the different methods, but there were still taxonomic similarities (excl. “EDTA”) (Figure 2b and Figure 3b) among the methods. This low level of overlap is probably due to the very low number of fungi in the samples. The most abundant fungal genera we found in the three methods “raw”, “pellet” and “total” were Aspergillus and Pseudeurotium. Aspergillus is a very diverse genus with over 250 species, which contains both pathogenic species, responsible for human and plant death, and species used in the food industry (Gibbons and Rokas 2013). However species of the genus Aspergillus are also known to contain collagenase producing fungi (Wanderly et al. 2017; Yakovleva, Khoang, and Nikitina 2006) which could be responsible for degradation of the collagen in the bone material. However, very little research has been published on the genetics of bone degrading fungi. The most abundant fungus genus found in the “EDTA” method were either yeast (cyberlindnera) (Kurtzman and Robnett 2013), pathogenic (Sarocladium) (Giraldo et al. 2015) or mushrooms (Clitocybe and Schizophyllum). This reflects the low resolution of the amplified ITS region and the fungal database used for taxonomic assignment. These genera that were present in the “EDTA” method were not present in any of the other methods. It could be speculated that the fungi present in the bone material will not be extracted when only using the supernatant as the rather large and heavy fungal cells will be spun down, and not incorporated in the further analysis. It could also be that the fungus we find in the supernatant is contamination in the EDTA solution and the reason we do not apply it in the other methods using EDTA is because it is in such low abundance it is overwhelmed by the bone fungus present in the pellet fraction.

Numerous studies have compared the efficacy of different DNA extraction kits in profiling microbial communities from different sample types (e.g. soil or saliva), and variations in the extracted communities are commonly reported between methods tested (Fredricks, Smith, and Meier 2005; Lazarevic et al. 2013; Leff et al. 1995; Santos et al. 2015; Vesty et al. 2017). Indeed, biases towards specific phyla are even reported when using different cell lysis methods (Lazarevic et al. 2013). In our data we observe variation within a single method (Figure 2a,b) and not just between methods. This could be due to the nonhomogeneous structure of the bone material itself, as other studies have shown a large variation in endogenous DNA yields between not only different skeletal parts (Mundorff and Davoren 2014), but also within the same bone (Pinhasi et al. 2015). It could also be due to the already mentioned low amounts of input DNA in the PCR reactions, which have been shown to have a significant impact on the sample profile (Kennedy et al. 2014). Contamination from the extraction kit, PCR reagents, laboratory equipment etc. can be an issue in low biomass microbial studies (Stinson, Keelan, and Payne 2018). However in our study both the extraction blank and the PCR blanks did not amplify, and the PCR blanks yielded no reads after sequencing. Although it would benefit the experiment to apply the different extraction approaches to a larger number of samples, our results still show that care should be taken when extracting microbial DNA from bones and the need to be consistent in whichever method chosen prior microbial DNA extraction. It also shows the importance of communicating which method that has been used in the study for other researchers to be able to compare results.

Conclusion

Overall, we found little overlap in the number of OTUs between the four different methods for both bacterial and fungal data. Fewer than 10% of the bacterial and 2% of the fungal OTUs were shared between the four different methods. The lack of shared OTUs among different molecular methods highlights the challenges and sensitivities of microbiome research.

Which method to choose depends on the research question asked. In a case where one wishes to obtain an overview of the microbial community present on a bone, disregarding a predigestion could be recommended. However to get full representation of the microbial community within the bone, whether the microbes associated with diagenesis or pathogens present in the animal’s blood system prior to death, then using the method where predigestion is used and both the supernatant and the pellet are extracted together could be recommended for both the fungal and bacterial communities. We do not recommend only using the supernatant after the predigestion step, prior to subsequent extraction of the bone material, as the fungal community showed very little overlap with the other methods as well as very little taxonomic diversity.

In this study, although the lowest amount of technical variability was observed in both the raw and total methods, the total method recovered microbes most relevant to diagenesis of mineralised tissue, and thus ultimately may prove to be the most suitable treatment of bone prior to DNA extraction for studies aiming to characterise the microbial community – something that we hope future studies may be able to explore in detail.

Acknowledgement

The authors would like to thank colleagues Sarah Mak, Kristine Bohmann, Katharina Fietz, Martin Nielsen, Christian Carøe and Shyam Gopalakrishnan from the EvoGenomics group at the GLOBE institute, University of Copenhagen, as well as two anonymous reviewers for very helpful guidance. We would also like to thank Conservator Kristian Murphy Gregersen from the Natural History Museum of Denmark for providing the cow humerus used in the experiment, and Henning Matthiesen, National Museum of Denmark for useful comments on the manuscript. This work was supported by the Danish Council for Independent Research under grant DFF –4090-00037.

Data availability statement

The data that support the findings of this study are openly available after publication of the paper.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Notes on contributors

Anne Marie Høier Eriksen is a Post Doc at the National Museum of Denmark and holds a PhD from the University of Copenhagen’s GLOBE Institute. Her research interests focus on the degradation processes associated with organic materials (e.g. bone and wood) in the marine environment, where she uses molecular approaches to studying the microbial community associated with diagenesis.

Lara Puetz is a PhD student at the University of Copenhagen’s GLOBE Institure who uses metagenomic tools to compare the microbiomes of different animal gut datasets in order to look at how community difference may affect phoenotypes of the host.

Carolina Rocha has a Bachelor in Agro-Genomic Sciences from The National Autonomous University of Mexico. Since she finished her bachelor thesis she has been working on microbiome projects, with a focus on microorganism especially pathogens. Currently she is studying pathogens from ancient individuals. In the future, she would like to study pathogens in modern individuals and study pathogen evolution.

Tue K. Nielsen is a post doc at the University of Copenhagen, Dep. of Plant and Environmental Sciences, Microbial Ecology and Biotechnology. He has experience with molecular microbial methods including high-throughput sequencing on 2nd- and 3rd-generation platforms, bioinformatics, and statistics. He has previously worked extensively with microbial community analysis using targeted amplicon sequencing, transcriptomics, and metagenomics. Furthermore, he has experience in working with comparative genomics of microbial genomes.

Lars H. Hansen is a professor at the University of Copenhagen, Dep. of Plant and Environmental Sciences, Microbial Ecology and Biotechnology. He uses molecular approaches to investigate complex questions in microbial ecology in both fundamental and applied sciences. His main fields of research include how bacteria perceive their environment, and change accordingly, via transfer and of bacterial plasmids as well as plasmid evolution in the context of bacterial resistance to antibiotics. In the past years, 2nd and 3rd generation DNA sequencing has become a favourite tool for investigating bacterial, bacteriophage and plasmid diversity in research related to the environment and bacteriophage-host interactions. His long-term experience in microbial sequencing projects has resulted in co-authoring more than 75 peer-reviewed publications on microbial sequencing in the last 5 years.

M. Thomas P. Gilbert is a professor at the University of Copenhagen’s GLOBE Institute. He principally applies ancient DNA methodologies to explore questions spanning evolution, ecology, archaeology and anthropology, but also in the context of degradation to both natural materials and the DNA that is within them. He additionally has an interest in environmental DNA studies, for use in both biomonitoring but also ecosystem reconstruction.

Additional information

Funding

This work was supported by the Danish Council for Independent Research under grant DFF –4090-00037