High-throughput sequencing to evaluate the effects of methamphetamine on the succession of the bacterial community to estimate the postmortem interval

Abstract In forensic medical examinations, estimating the postmortem interval (PMI) is an important factor. Methamphetamine (MA) is a synthetic stimulant that is commonly abused, and estimation of the PMI after MA abuse has become one of the main tasks in forensic investigation. Microorganisms play a vital role in carrion decomposition. Analysing the bacterial succession patterns can be used as a forensic tool to estimate the PMI. The present study aimed to analyse bacterial succession changes during the decomposition of MA to estimate the PMI. We analysed bacterial communities in rabbits treated with three different concentrations of MA (0, 22.5, and 90 mg/kg) under the natural conditions of 20 °C and 70% humidity by sequencing 16S rRNA gene amplicons using the Illumina MiSeq system. We obtained 2 374 209 high-quality sequences and 2 937 operational taxonomic units (OTUs). The relative abundances of the bacterial communities varied markedly in response to different MA concentrations. Interestingly, in response to the different concentrations of MA, Bacteroidetes became disparate in the rectum in the late PMI. Increased numbers of bacterial taxa were identified in the rectum and buccal cavity samples, except at the highest concentration of MA in the rectum samples when PMI was 0–h, than were present in live rabbits. Meanwhile, the PMI correlated significantly with bacterial succession at different taxonomic levels. Our results suggested that bacterial community succession could be used as a “microbial clock” to estimate the PMI in cases of MA-related death; however, further study is required to gain a deeper understanding of this concept.


Introduction
Estimation of the postmortem interval (PMI) is an important and difficult aspect of forensic investigation. Generally, the PMI refers to the period from when death occurs to when the corpse is found, which can also be called the "postmortem experience time" [1,2]. Currently, most methods to estimate PMI rely on the morphological method, including early postmortem phenomena (algor mortis, rigor mortis, and livor mortis) and late postmortem phenomena (greenish discolouration, putrefactive blister, and putrefactive networks) [3][4][5]. However, decomposition is a complex process, which is greatly affected by biological (cellular enzymes, microbes, and insects) and environmental (temperature, weather, and humidity) factors [6], and can only produce a rough estimate of the PMI. In addition, other methods have been used to estimate the PMI, such as the developmental or evolutionary patterns of necrophilic insects [7][8][9][10][11][12] and genomic and transcriptomic methods [13][14][15][16]. However, necrophilic insect activities are affected by sunlight, temperature, and weather, which might lead to inaccurate estimates of the PMI [17]. Similarly, genomics and transcriptomics are affected by multiple circumstantial and environmental factors, in which the accuracy and precision decrease as the PMI increases. Therefore, despite decades of research, the accuracy of estimating the PMI has not improved significantly, and no method can be used reliably to estimate the PMI accurately [18]. To overcome the limitations in practice, the great potential of the microbiome as a tool for estimating the PMI in forensic has emerged.
Microorganisms are a form of ubiquitous evidence that present regardless of the season, which often responds to changes in their environment in a predictable way. The postmortem microbial community's response to PMI is predicable [19][20][21]. Specifically, Metcalf et al. [22] reported the use of a mouse model that allowed them to propose the concept of a "microbial clock", indicating that the microbial community decomposition of cadavers might be developed as a forensic tool to assess the PMI. Similarly, Pechal et al. [20] used a swine carcasses model to reveal a significant negative linear relationship between the microbial phyla and the richness of family taxa during the decomposition progresses. In summary, these studies indicated that the succession of microorganisms during the decomposition processes seems to be a foreseeable process, which has a significant impact on estimating the PMI in forensic practice. Previous studies have shown that during PMI estimation, drugs and poisons could change the direction of microorganisms succession [23][24][25]. Therefore, when estimating the PMI via changes in the microbiome, it is necessary to consider the effects of drugs or poisons on the succession of microbial communities.
Methamphetamine (MA) is a synthetic stimulant that is commonly abused in many countries worldwide. Generally, the phenomenon of MA abuse has become a global problem. With the prevalence of MA abuse, negative socioeconomic, cognitive, and legal consequences have followed [26]. Therefore, in cases of suicides or accidents, estimation of the PMI of MA abuse has become one of the main tasks in forensic investigation. In recent years, studies have shown that MA could change the development time of the larvae of necrophilic insects and cause deviations in the estimated PMI [5,27,28]. Among these, MA can change the balance of intestinal microbes. Ning et al. [29] proved that the use of MA caused changes in the intestinal flora of rats in a study of MA-induced conditioned place preference rats. Chen et al. [30] revealed that MA promotes intestinal inflammation by increasing the relative abundance of pathogenic bacteria in the intestine and reducing the expression of tight junction (TJ) proteins in the intestines. Although studies have shown that the abuse of MA can lead to significant time-and structure-dependent changes in the diversity of the gut microbiome and taxonomic structure, only a handful of studies have described the relationship between the effect of MA on the microbiome and PMI estimation.
In the present study, we aimed to explore the relationship between the changes in the buccal cavity and rectum microbiome for estimating the PMI under three different concentrations of MA using high-throughput sequencing. The results revealed that at different taxonomic levels, the succession of microorganisms is related significantly to the PMI for estimating the use of MA.

Study sites and grouping of samples
As models of human decomposition, we used adult female New Zealand rabbit carcasses to study the changes in the bacterial community during decomposition under three different concentrations of MA. The experiments were conducted in December, 2019 in Changsha City, Hunan Province, Central South China. The laboratory temperature was approximately 20 °C with 70% humidity and natural light. The rabbits, each weighing 2.0 kg, were divided randomly into three groups. In the control group (Group A), 0.6 mL saline solution was injected through the ear margin veins, and death was caused by air embolization through the ear margin veins at 1.5 h after injection. Two different doses (22.5 and 90 mg/kg) of MA, diluted in 0.6 mL of saline solution, were injected into the ear margin veins. MA at 22.5 mg/kg represented the half lethal concentration (Group B), and in this group, death was caused by air embolization through the ear margin veins at 1.5 h after injection. MA at 90 mg/kg represented twice the lethal concentration (Group C) and caused immediate death. Each corpse was individually put into plastic boxes and placed into the centre of the laboratory. The boxes contained sliver sand to absorb the putrefactive liquid to avoid contamination of the surroundings. To simulate human decomposition, the necrophagous insects had access to the corpses in the natural environment. Each day, the physical changes in the carcasses were observed and the surface changes were recorded.

Collecting samples of the bacterial community
We sampled the bacterial communities before death (5 min before killing), immediately after death (less than 10 min postmortem), and at 2 h, 6 h, 12 h, 24 h, 3 d, 5 d, 10 d, and 20 d postmortem. Sterile cotton applicators scrubbed with sterile ultrapure water were used to swab the rectum and buccal cavity lightly for 1 min. The cotton applicator was transferred to a 1.5 mL sterile Eppendorf tube with 1 mL sterile ultrapure water. Samples were then snap frozen at −80 °C for later use.

Extraction of DNA, PCR amplification, and sequencing
Total bacterial DNA was extracted from the New Zealand rabbit buccal cavity and rectum samples using a MoBio PowerSoil DNA Isolation Kit (Mo Bio Laboratories, Carlsbad, CA, USA), following the manufacturer's instructions. The concentration and purity of the DNA were assessed using 1% agarose gel electrophoresis. In 1.5 mL sterile Eppendorf tubes, the DNA samples were diluted using sterile water to 1 ng/mL. In each group, the DNA samples from the rabbit carcasses were mixed in equal concentrations at each time point. Novogene Biological Information Technology Co. (Beijing, China) completed the sequencing operation using Illumina MiSeq sequencing.
A 515F/806R primer set 50 was used to conduct the PCR reactions [31], amplifying the 16S rRNA gene V3-V4 hypervariable regions (forward primer: 5′-GTGCCAGCMGCCGCGGTAA-3′, reverse primer: 5′-GGACTACHVGGGTWTCTAAT-3′). The reverse primer included a 6-bp error-correcting barcode that was unique for each sample. The PCR reactions comprised 15 μL of Phusion High-Fidelity PCR Master Mix (New England Biolabs, Ipswich, MA, USA), 0.2 μmol/L forward and reverse primers, and 10 ng of template DNA in a 30-μL reaction volume. The PCR reactions were carried out as follows: 98 °C for 1 min; 30 cycles of 98 °C for 10 s, 50 °C for 30 s, and 72 °C for 30 s; a final 5 min extension at 72 °C. Agarose gel electrophoresis (2%) was used to confirm the resulting amplicons. A GeneJET Gel Extraction Kit (Thermo Fisher Scientific, Carlsbad, CA, USA) was used to purify the 400-450 bp amplicons. The amplicons were quantified, and equal amounts of the purified amplicons were combined into a single tube. An NEB Next Ultra DNA Library Prep Kit for Illumina (New England Biolabs) was used to generate the sequencing libraries, according to the suppliers protocols, including the addition of index codes. A Qubit 2.0 Fluorometer (Thermo Fisher Scientific) and an Agilent Bioanalyzer 2100 system (Agilent, Santa Clara, CA, USA) were used to assess the quality of the libraries. The Illumina MiSeq platform (Illumina, San Diego, CA, USA) was used for sequencing, generating 300 bp paired-end reads.
FLASH version 1.2.7 was used to merge the paired-end reads from the original DNA fragments [32]. Low-quality parts were removed using QIIME version 1.7.0, and the primer sequences and barcodes were trimmed produce the original reads [33]. QIIME was used to filter out sequences ambiguous bases (N) or low-quality bases [34] and UCHIME was used to remove chimeric sequences to obtain clean reads [35]. The unique barcodes on the remaining sequences were used to assign them to the samples, and before statistical analysis, the sequences were subjected to rarefaction at the level of 9200.
The UPARSE pipeline (version 7.0.1001) was used to cluster the sequences [36], and a threshold of 97% similarity was used to assign similar sequences to operational taxonomic units (OTUs). Among the similar sequences, a representative sequence from each OTU was identified, which was the longest sequence with the maximum number of hits to the other representative sequences. Using a minimum identity of 80%, the representative sequences were aligned at the GreenGenes database [37]. PyNAST software (version 1.2) was used to determine the phylogenetic relationships among the representative sequences, and for rapid multiple sequence alignment [31], and the GreenGenes database "Core Set" data [38] were used. "Rare taxa" were defined as those taxa with a relative abundance <1% of the total abundance of all samples. Unclassified reads were categorised as "others".
Three phylogenetic diversity metrics were used to assess the α-diversity: the Shannon index, Chao 1, and Observed species. According to Observed species, we generated rarefaction curves. Both unweighted and weighted UniFrac distances were used to evaluate the β-diversity between the bacterial communities [39]. The unweighted pair group method with arithmetic mean (UPGMA) was used to stratify the hierarchical clustering of the samples. The UniFrac distances of weighted and unweighted distance matrices were subjected to principal components analysis (PCoA) to visualise the difference in bacterial community structure and composition, respectively. Based on the Bray-Curtis dissimilarity distance matrices, analysis of similarities description was performed to examine the differences in the community composition between the groups of samples [40].

Results
The study was performed at a daily temperature of approximately 20 °C with 70% humidity, in a reversed light cycle of 12 h:12 h light:dark. The decomposition process of the control and treatment groups showed no obvious differences according to the physical changes of the carcasses. In the current study, we observed five stages of decomposition according to the physical changes: fresh, bloated, active decay, advanced decay, and dry remains. In the first stage of fresh decomposition (0-1 d), the necrophagous insects arrived at the carrion immediately, indicating that the presence of MA in the rabbits did not impede their activities. Within the bloated stage of death (2-3 d), necrophagous insect eggs were oviposited on eyes and noses of the cadavers. During the active decomposition stage (4-5 d), a strong odour emerged and putrefactive fluid leakage was observed. The characteristics of the advanced decomposition stage (6-10 d) were the lack of odour, the removal of most of the soft tissues, and a high number of necrophagous insect larvae. Decomposition finished with the dry remains (11-20 d) and during this stage, we observed post-feeding larvae and pupa of the necrophagous insects.
Thirty buccal cavity samples and 30 rectum samples were collected at three different MA concentrations (0, 22.5, and 90 mg/kg). Sequencing was successful for all sample. After quality filtering, failed sample and low numbers of sequences removal, 2 545 378 raw sequences were generated, among which we retained 2 374 209 high-quality sequences with an average length of 252 bases. For each sample, the mean number of sequences was 4 242. The species representation in each sample had approached the plateau phase according to the rarefaction curves, thus additional sequencing effort was unlikely to detect more bacteria (Figure 1). The UPARSE pipeline clustered the high-quality sequences into 2 937 OTUs at a 97% identity threshold.
Principal coordinate analysis (PCoA) of the clustering patterns of samples from the rectum and buccal cavity during decomposition provided the decomposition pattern in a two-dimensional space for subsequent weighted UniFrac distance analysis ( Figure 3). The weighted UniFrac PCoA analysis showed that rectal samples collected from rabbits up to 1 d formed a unique cluster that was separated from the samples from the buccal cavity. Principal coordinate 1 (PC1) and Principal coordinate 2 (PC2) analysis (explaining 20.08% and 14.34% of the variance, respectively) indicated that the rectal and buccal microbial communities were separated during decomposition. The hierarchical clustering of buccal cavity and rectum samples using UPGMA permitted analysis of the samples according to their weighted UniFrac matrix. Rectal samples collected from rabbits until 1 d formed a unique cluster that was separated from the samples from the buccal cavity. This result was similar to that of the weighted UniFrac PCoA. Previous studies have demonstrated that the buccal cavity community differs from that of the rectum of cadavers during the expansion phase of Figure 1. rarefaction curves of observed species number clustered at 97% sequence identity across all samples. sample named to refer to samples as described in Table 1.
decay [41]. In the present study, the buccal cavity and rectal microbial communities were distinguished during decomposition, suggesting that the sequencing data were reliable (Figure 4).
The 16S RNA sequences were then classified at the phylum and the family levels to identify the succession of the bacterial community structure during decomposition. The relative abundance of different microflora showed obvious changes during the decomposition of the buccal cavity and rectum of rabbit carcasses (Figures 5 and 6). In the buccal cavity ( Figure 5A-C), the dominant phylum was the Proteobacteria (63.80% ± 18.32% (SD) in Group A, 65.44% ± 14.92% in Group B, and 62.38 % ± 29.04% in Group C) with the progression of decomposition, except in the immediately after death sample of Group A, in which the Bacteroidetes were the dominant taxon (55.10%). Similarly, there was a large difference between the buccal cavity samples of the experimental group and the control group. At 6 h  35.10%, and 1.01%, respectively) and Group B (0.56%, 1.34%, and 0.01%, respectively), while they showed a tendency to increase in Group C (0.07%, 0.37%, and 1.06%, respectively). The rabbit's corpse is an energy body for microorganisms, however, as corruption progresses, the body's energy gradually decreases. These data showed that Bacteroidetes were more affected by MA, and the higher the concentration, the stronger the inhibitory effect on the Bacteroidetes. Interestingly, the numbers of Actinobacteria decreased gradually with the progression of decomposition, and almost vanished from 5 d postmortem onwards in all groups. At the family level in the buccal cavity ( Figure 5D-F), as decomposition progressed, the numbers of Xanthomonadaceae increased gradually, becoming the dominant taxon at 10 d postmortem in Group B (63.02%) and  ) and Group c (orange dot) were represented. sample names refer to samples as described in Table 1. communities. sample names refer to samples as described in Table 1. The difference between the samples treated with MA at 22.5 mg/kg and 90 mg/kg were tested using analysis of similarities description. Between the two groups, we noted differences for several taxa. In the buccal cavity, Capnocytophaga until 1 d postmortem and Elizabethkingia in the late PMI had higher prevalences in Group B than in Group C. Indeed, Capnocytophaga and Elizabethkingia were not detected in Group c ( Supplementary Figure 2A-C). In the rectum, Carnobacterium had a higher prevalence in Group b than in Group c, and Anaerococcus was only detected in Group c (Supplementary Figure 2D-F).

Discussion
In this study of the postmortem microbiome in response to MA, in the buccal cavity of live rabbits, the dominant phyla were Proteobacteria, Bacteroidetes, Firmicutes, and Actinobacteria, and the dominant genera were Actinobacillus and Riemerella. In the rectum of live rabbits, the dominant phyla were sample names refer to samples as described in Table 1.
Bacteroidetes, Firmicutes, and Proteobacteria, and the dominant genera were Bacteroides and Corynebacterium. We chose the buccal cavity and rectum as the sampling point for microbes because of the convenience of the sampling process, which simplified subsequent analysis. The dominant phyla and genera in the live rabbits were similar to those in samples from the mouth and stool of healthy humans, according to the Human Microbiome Project (HMP) [42]. The HMP reported that the microbiome of the human stool is rich, whereas that of the oral cavity is more limited [42][43][44], which was similar to our results from samples collected from live rabbits.
New Zealand rabbits not only have similar decay processes to humans, but also harbour similar microbial communities. Among this, similar to research using rat models, the use of rabbit models permitted replicate experiments incorporating a large sample size to be performed. This not only minimized any possible experimental error, but also allowed the assessment of the extent of intra-individual microbiota variation in microbiota that occurred during decomposition [3,6], which indicating that these rabbits are suitable models to study human decomposition.
Compared with that from live rabbits, more bacterial taxa were identified in the buccal cavity immediately after death. The mouth contacts the outside environment directly, and thus external factors, such as humidity and temperature, can easily effect changes in the microbiome. Therefore, to reduce the effects of the external environment to a minimum, we controlled the temperature, humidity, and other controllable factors. In particular, in the rectal samples, there were marked differences pre-and immediately postmortem. Previous studies found no obvious differences in samples collected pre-and immediately postmortem from Sprague Dawley rats  Table 1. [3]. This apparent discrepancy might be caused by species differences between the experimental animals; however, the effect of different concentrations of MA on the experimental animals cannot be ruled out.
According to previous studies, in the bloat stage of decomposition and in healthy humans, the detected oral microbial communities were derived from the gastrointestinal tract communities [41,42]. However, interestingly, we noted that samples from the buccal cavity and rectum were separated in the early PMI, but were clustered together in the late stage. Previous experimental data revealed that after death, saliva, which generates antibacterial activity, ceases to be produced and the intestinal microenvironment can change, thus the limiting factors of microorganisms are gradually weakened [45,46]. Furthermore, with increasing PMI, the differences between the buccal cavity and rectum microorganisms gradually decreased. During the decomposition process, the abundance of bacteria increased, but their richness decreased, which could be used for PMI estimation. In addition, both the rectal and buccal cavity samples contained reads from unclassified organisms that were within acceptable limits. Previous studies revealed that in metagenomic sequencing, approximately 80% of the identified bacteria communities were nonculturable [44]. Although some of the unclassified sequences might represent PCR errors or sequencing artifacts, the large numbers of unclassified sequences indicated the presence of novel bacterial taxa. The presence of novel species can be inferred when the phenotypical characteristics and/or 16S rRNA sequences of unknown bacteria (greater than 0.5% difference) Figure 6. Bacterial community structure variation during decomposition in the rectum at the phylum level and the family level. relative abundance of bacterial phyla during decomposition in the rectum in Group a (a), rectum in Group b (B), and rectum in Group c (c). relative abundance of bacterial families during decomposition in the rectum in the Group a (D), the buccal cavity in Group b (e), and the buccal cavity in Group c (F). sample names refer to samples as described in Table 1.
are significantly different from those of the most closely related bacteria. Besides, the bacteria detected on and in a cadaver can be influenced by many factors. For example, the individual's initial microbiome, the environment surrounding the decomposing cadaver, the method of sample collection, and even variations in swab pressure and the number of swabs used during sample collection.
Our findings revealed that during decomposition, the relative numbers of dominant phyla, i.e. Proteobacteria, Bacteroidetes, Firmicutes, and Actinobacteria, varied significantly, which might assist the estimation of the PMI in cases of MA abuse. At 3 d after death, the Bacteroidetes were prevalent members of the bacterial community in the buccal cavity, becoming even more abundant compared with their level in the control groups (Group A and Group B). Interestingly, the most prevalent phylum on day 10 and day 20 postmortem was the Proteobacteria, which becoming more abundant compared with their levels in the control groups (Group A and Group B). Meat spoilage is commonly associated with the Proteobacteria, and these bacteria have been detected in the skin of slaughtered animals [47]. MA exposure increased the content of conditionally pathogenic bacteria with pro-inflammatory effects, such as the Proteobacteria [30,[48][49][50]. Analysis of the groups with different concentrations of MA revealed that the abundance of the Firmicutes increased during decomposition. The two main phyla in the human gut are reported to be the Bacteroidetes and Firmicutes. In addition, our determination of the relationship between MA administration and intestinal flora revealed that the Firmicutes accounted for the majority of the changed taxa and the propionate-producing genus Phascolarctobacterium was repressed by MA administration [29]. These findings indicated that MA was closely related to changes in the intestinal microbiome. Actinobacteria are widely distributed in terrestrial and aquatic ecosystems, especially in soil, where they play a significant role in the recycling, decomposition, and the formation of refractory biomaterial [51]. Our results further indicated that the Gammaproteobacteria class was predominant in all groups in the rectum and buccal cavity during the late stage of decomposition. Gammaproteobacteria have different metabolic capabilities and can breakdown more complex molecules [52]. Our findings supported the hypothesis that the notable variation of microbial communities might aid PMI estimation in the presence of MA abuse.

Conclusion
In summary, we sought to explore the relationship between the buccal cavity and rectum microbial changes and the PMI under three different concentrations of MA using high-throughput sequencing technology.
Our results showed that MA induced significant changes in bacterial community succession by decreasing the abundance of probiotics and increasing the abundance of conditioned pathogens with pro-inflammatory effects. As expected, the results support the theory that significant changes in the microbial community might contribute to the estimation of the PMI in the context of MA abuse. Taken together, although further investigation is needed, these findings provide new ideas for the estimation of the PMI in the context of MA abuse using changes in the microbiome. Our research had some limitations. Particularly, despite verifying the effect of MA on the microbial community, we have not thoroughly studied the underlying mechanism of the link between the change of the microbial community caused by MA and the estimation of the PMI. Further experiments are expected to clarify the potential connection between MA-induced microbial community changes and PMI estimation.