Postmortem submersion interval prediction model based on the rat muscle microbiome

TL;DR

This study explores how the microbiome in rat skeletal muscle changes predictably after submersion in water, offering a new method to estimate post-mortem submersion intervals.

Contribution

The study introduces a high-precision two-stage machine learning model for PMSI estimation using skeletal muscle microbiome data, independent of cause of death.

Findings

The rat skeletal muscle microbiome shows a predictable successional pattern after submersion.

A two-stage model achieved 90.9% accuracy in classifying early and late phases of submersion.

The model's mean absolute error was 0.303 days in the early phase and 1.293 days in the late phase.

Abstract

Accurate estimation of post-mortem submersion interval (PMSI) is a critical challenge in forensic science. The current research has largely focused on the microbial communities in the skin and gut, which are susceptible to environmental contamination, while the potential of internal tissues remains underexplored. This study aimed to investigate whether the microbiome of skeletal muscle, a relatively closed ecosystem, undergoes a predictable succession following submersion in water PMSI and to evaluate its potential for building a high-precision PMSI prediction model, which is independent of the cause of death (drowning vs. post-mortem submersion).

Using 72 male Sprague–Dawley rats, we established drowning (D group) and post-mortem submersion (PS group) models. After submersion in natural aquatic environment for 14 days, skeletal muscle samples were collected at six time. The microbial communities were profiled by high-throughput sequencing of the V3–V4 region of the 16S rRNA gene, followed by analyses of alpha and beta diversity. Based on the observed successional patterns, a two-stage prediction model combining classification and regression algorithms (e.g., random forest, RF) was developed.

The skeletal muscle microbiome exhibited a significant and predictable successional pattern, clearly partitioning into an early-phase (0–3 days) and a late-phase (5–14 days) (PERMANOVA, p < 0.001). This succession was characterized by a shift from the dominant community of Proteobacteria to the dominant community of Firmicutes. Importantly, the cause of death did not significantly impact either the alpha or beta diversity of the microbial communities (PERMANOVA, p = 0.251). The resulting two-stage prediction model demonstrated excellent performance: the classification model distinguished the early and late phases with an accuracy of 90.9% (AUC = 0.9504), and the mean absolute errors (MAE) of regression models was 0.303 days in the early phase and 1.293 days in the late phase.

The rat skeletal muscle microbiome undergoes a regular and predictable post-mortem succession unrelated to the cause of death. The stable “microbial clock” within the internal tissue allows the construction of a high-precision two-stage machine learning model for PMSI estimation. Our results establish skeletal muscle as a highly promising new target for forensic microbiology, offering a robust theoretical basis and technical approach to resolving challenges in long-term PMSI estimation.

Figures

Compositional dynamics and successional patterns of the muscle microbial community at different taxonomic levels. (A) Stacked bar chart showing the relative abundance of dominant bacterial phyla over the 14-day PMSI. (B) Stacked bar chart displaying the relative abundance of dominant bacterial genera. The legend shows the top 10 most abundant genera. (C) Heatmap illustrating the Z-score normalized relative abundance of the top 30 most abundant genera across all samples. Samples are ordered by PMI in days, and genera are clustered hierarchically based on their abundance profiles. (D) Sankey diagram illustrating the compositional flow from sample groups (ordered by PMSI) to phylum and genus levels. The width of the flows is proportional to the relative abundance.

Comparison of muscle microbial community diversity between the drowning (D) and post-mortem submersion (PS) groups. (A) Principal co-ordinates analysis (PCoA) plot based on Bray–Curtis dissimilarity, illustrating the overall microbial community structure. The first two principal coordinates explain 22.82 and 13.91% of the total variation, respectively. Each point represents an individual sample, colored by group (D: red; PS: blue). Ellipses indicate the 95% confidence interval for each group. (B,C) Boxplots comparing the alpha diversity between the D and PS groups at each PMSI. Community richness was estimated using the Chao1 index (B), and community diversity was estimated using the Shannon index (C). No significant differences (ns, p > 0.05) were observed between the two groups at any time point according to the Wilcoxon rank-sum test. The boxes represent the interquartile range (IQR), the central line indicates the median, and the whiskers extend to 1.5 times the IQR.

Distinct microbial community structures in rat muscle at early and late post-mortem submersion intervals. (A) Principal co-ordinates analysis (PCoA) and (B) non-metric multidimensional scaling (NMDS) plots illustrating the beta-diversity of microbial communities, based on the Bray–Curtis dissimilarity metric. The first two principal coordinates in the PCoA plot explain 22.8 and 13.94% of the total variation, respectively. The NMDS plot has a stress value of 0.1708. Each point represents an individual sample. Samples are colored according to their post-mortem submersion interval (PMSI) phase: yellow for the early-phase (EP; 0–3 days) and purple for the late-phase (LP; 5–14 days). Ellipses represent the 95% confidence interval for each group. The separation between the two phases was statistically significant (PERMANOVA, p < 0.001).

Performance evaluation of machine learning models for classifying early-phase (EP) and late-phase (LP) post-mortem submersion intervals. (A) Confusion matrices for four different classification models: random forest (RF), logistic regression (LR), k-nearest neighbors (KNN), and support vector machine (SVM). The values in each cell represent the number of samples in the test set. The color intensity corresponds to the number of samples. (B) Receiver operating characteristic (ROC) curves for the four models. The area under the curve (AUC) for each model is displayed in the legend.

Performance of the two-stage random forest regression models for PMSI prediction. Scatter plots showing the correlation between the observed PMSI (in days) and the PMSI predicted by the regression models for (A) the early-phase (EP; 0–3 days) and (B) the late-phase (LP; 5–14 days). Each point represents an individual sample, colored by its inclusion in either the training set (blue) or the test set (red). The solid black line represents the ideal prediction line (y = x). The mean absolute error (MAE) for both the training and test sets is displayed in the inset box for each model.

Full text

Introduction

1

The accurate estimation of the post-mortem submersion interval (PMSI) is a cornerstone of forensic science practice, providing critical information for criminal investigations, suspect validation, and event timeline reconstruction (van Daalen et al., 2017; García et al., 2021). Traditional methods, such as assessing physical changes like body temperature, rigor mortis, and livor mortis, are primarily applicable to the early post-mortem period (typically within 72 h) and their accuracy is easily compromised by various factors including ambient temperature and individual variations (Pittner et al., 2020; Zissler et al., 2020). Furthermore, forensic entomology provides limited information for PMSI estimation of submerged remains, as colonization by terrestrial insects typically only occurs once the cadaver floats to the surface. Consequently, the development of objective, universally applicable methods for estimating longer-term PMSI remains a significant challenge and an urgent need in the forensic field.

In recent years, the advent of high-throughput sequencing technologies has catalyzed the emergence of forensic microbiology, opening a novel frontier for PMSI estimation. A substantial body of research has demonstrated that upon the host’s death, its symbiotic microbial communities, both internal and external, undergo dynamic and predictable successional changes (Metcalf et al., 2013; Metcalf et al., 2016). This concept of a “microbial clock” has been validated in various animal models and human cadaver studies. Current research has predominantly focused on communities in the gut (Zhang et al., 2022; DeBruyn and Hauther, 2017), oral cavity (Metcalf et al., 2017), and on the skin (Lax et al., 2015), collectively establishing a strong time-dependent relationship between microbial succession and PMSI and enabling the development of preliminary prediction models using machine learning algorithms (Liu et al., 2020).

However, despite this significant progress, two core challenges hinder the translation of these findings into routine forensic practice. First, the focus on the gut and body surfaces—sites in direct contact with the external environment—means that these microbial communities are highly susceptible to random invasion by environmental microbes, contamination from the surrounding medium (e.g., water, soil), and disturbances from insect activity, all of which can compromise the stability and accuracy of the “microbial clock” (Bone et al., 2024; Wallace et al., 2021). Second, the specific impact of the cause of death on microbial succession remains largely unclear. Different manners of death, such as acute trauma versus drowning, may result in vastly different initial physiological states (e.g., circulatory status, tissue oxygenation), but whether this systematically alters the subsequent trajectory of microbial succession remains an unresolved question.

To address these knowledge gaps, this study focuses on a long-overlooked yet unique target for forensic microbiology: skeletal muscle. As an internal tissue with immense biomass that is relatively enclosed during the early stages of decomposition, skeletal muscle may provide a more stable source of information for PMSI estimation, one less affected by random external factors. We, therefore, hypothesized that the microbial community within skeletal muscle undergoes a regular and predictable succession following submersion in water (PMSI) and that this internal successional pattern remains highly consistent regardless of the cause of death (drowning vs. post-mortem submersion). In this study, we employed a controlled rat model to systematically track the dynamic changes in the skeletal muscle microbiome over a 14-day submersion period. By integrating machine learning algorithms, we aimed to answer the following questions: (1) Does a clear, time-dependent successional pattern exist in the skeletal muscle microbiome? (2) Does the manner of entry into the water significantly affect this pattern? (3) Can a PMSI prediction model based on skeletal muscle microbial data achieve the accuracy required for forensic applications? Our results reveal a distinct succession characterized by a shift from Proteobacteria to Firmicutes dominance and demonstrate the successful construction of a high-precision, two-stage PMSI prediction model that is robust to the cause of death, offering a novel and potent solution for forensic PMSI estimation.

Materials and methods

2

Animal model and experimental design

2.1

All animal experiments were conducted in strict accordance with animal welfare and ethical guidelines and were reviewed and approved by the Institutional Animal Care and Use Committee of Guangzhou Hua Teng Biotechnology Co., Ltd. (Approval No. B202503-1). This study utilized 72 healthy adult male Sprague–Dawley (SD) rats, aged 8 weeks, sourced from Zhuhai BesTest Bio-Tech Co., Ltd. All rats were acclimatized for 1 week in a barrier facility (temperature 20–26 °C, relative humidity 40–70%, 12 h/12 h light/dark cycle, pressure differential ≥10 Pa, air changes ≥15 per hour), with ad libitum access to irradiated standard chow (Jiangsu Co-Shine Bio-Tech Co., Ltd.) and autoclaved drinking water.

To investigate the influence of the cause of death on post-mortem microbial succession, the 72 rats were randomly assigned to either a drowning group (D, n = 36) or a post-mortem submersion group (PS, n = 36). The experiment was conducted under natural outdoor conditions with a recorded ambient temperature ranging from 24.8 °C to 29.6 °C. To ensure environmental consistency, the D group and PS group rats were housed in two separate large transparent water tanks, filled with natural water to a depth of over 30 cm. For the D group, rats were placed in mesh bags and subjected to repeated cycles of submersion and air exposure to simulate drowning until vital signs ceased. For the PS group, rats were euthanized via CO₂ inhalation. Following confirmation of death, carcasses from both groups, still within their mesh bags, were fully submerged in the water tanks, with each carcass housed individually. To observe the dynamic changes, both the D and PS groups were further divided into six temporal subgroups corresponding to post-mortem submersion intervals of 0, 1, 3, 5, 7, and 14 days, with each time point containing six biological replicates.

Sample collection and preservation

2.2

At each designated time point, all six rats from the corresponding subgroup were retrieved for sample collection. All procedures were performed under sterile conditions, using a separate set of autoclaved dissection instruments for each animal. The skin and fascia of the hind limb were aseptically incised to expose the quadriceps or gastrocnemius muscle. A tissue block of approximately 0.5–1.0 g was excised from the deep central region of the muscle. All collected samples were immediately placed in pre-labeled, cryo-resistant screw-cap tubes, flash-frozen in liquid nitrogen, and subsequently stored at −80 °C until DNA extraction.

Microbiome DNA extraction and 16S rDNA sequencing

2.3

Total microbial DNA was extracted from approximately 0.5–1.0 mg of each skeletal muscle samples using the Bacterial Genomic DNA Extraction Kit (DP302-02, Tiangen Biotech (Beijing) Co., Ltd.), strictly following the manufacturer’s protocol. The concentration and purity of the extracted DNA were quantified using a Qubit^®^ 4.0 Fluorometer (Invitrogen, United States) with the dsDNA HS Assay Kit and assessed for integrity via 1% agarose gel electrophoresis.

The V3–V4 hypervariable region of the bacterial 16S rRNA gene was amplified using the primers 341F (5′-CCTACGGGNGGCWGCAG-3′) and 805R (5′-GACTACHVGGGTATCTAATCC-3′), synthesized by Sangon Biotech (Shanghai) Co., Ltd. PCR was performed in a 25 μL reaction volume containing 12.5 μL of Phusion Hot start flex 2X Master Mix, 2.5 μL of each 10 μM primer, and 50 ng of template DNA. The thermal cycling conditions were: an initial denaturation at 98 °C for 30 s; followed by 32 cycles of 98 °C for 10 s, 54 °C for 30 s, and 72 °C for 45 s; and a final extension at 72 °C for 10 min.

Amplified PCR products were verified by 2% agarose gel electrophoresis and purified using AMPure XT beads (Beckman Coulter Genomics, United States). The quality of the resulting libraries was assessed with an Agilent 2100 Bioanalyzer (Agilent, United States). Qualified libraries were pooled in equimolar amounts and sequenced on an Illumina NovaSeq 6000 platform using a NovaSeq 6000 SP Reagent Kit (500 cycles) for PE250 paired-end sequencing. This sequencing service was provided by LC-Bio Technology Co., Ltd. (Hangzhou, China).

Bioinformatic and statistical analysis

2.4

Raw sequencing data were processed using the QIIME2 pipeline (v2019.7). After demultiplexing, primers and barcodes were removed using the cutadapt plugin. The DADA2 plugin was then employed for sequence denoising, merging of paired-end reads, chimera removal, and generation of an amplicon sequence variant (ASV) feature table. The amplified V3–V4 region resulted in merged sequences of approximately 460 bp. A total of 7,587,720 high-quality reads were obtained from all samples after filtering. Rarefaction curves confirmed that the sequencing depth was sufficient to capture the majority of microbial diversity in each sample (Supplementary Figure 1). Taxonomic classification of ASVs was performed against the SILVA database (Release 138).

All statistical analyses were conducted in R (v4.1.3). Alpha diversity indices (Chao1, Shannon) and beta diversity metrics (Bray–Curtis) were calculated using the vegan and phyloseq packages. The Wilcoxon rank-sum test was used to compare alpha diversity between groups. Permutational multivariate analysis of variance (PERMANOVA) was employed to test for significant differences in community structure. Machine learning models were constructed using the randomForest package, with the relative abundances of microbial genera as predictor variables and the PMSI as the response variable. The performance of classification models was evaluated by the area under the receiver operating characteristic (ROC) curve (AUC), while regression model performance was assessed by the mean absolute error (MAE). A p-value <0.05 was considered statistically significant for all tests.

Results

3

Dynamic changes in the taxonomic composition of the microbial community

3.1

To identify the key microbial taxa driving community succession, we analyzed the relative abundance of species at different taxonomic levels. At the phylum level, the microbial community composition underwent a systematic change with increasing post-mortem submersion interval (PMSI). At the initial time point (day 0), the community was primarily composed of Firmicutes and Proteobacteria. As time progressed, the relative abundance of Proteobacteria consistently increased, rising from less than 20% initially to over 60% by day 14, while the relative abundance of Firmicutes steadily decreased (Figure 1A).

Compositional dynamics and successional patterns of the muscle microbial community at different taxonomic levels. (A) Stacked bar chart showing the relative abundance of dominant bacterial phyla over the 14-day PMSI. (B) Stacked bar chart displaying the relative abundance of dominant bacterial genera. The legend shows the top 10 most abundant genera. (C) Heatmap illustrating the Z-score normalized relative abundance of the top 30 most abundant genera across all samples. Samples are ordered by PMI in days, and genera are clustered hierarchically based on their abundance profiles. (D) Sankey diagram illustrating the compositional flow from sample groups (ordered by PMSI) to phylum and genus levels. The width of the flows is proportional to the relative abundance.

At the genus level, we observed more detailed successional patterns. A heatmap analysis revealed two main successional clusters based on the hierarchical clustering of genus abundance patterns (Figure 1C). During the early decomposition phase (0–3 days), genera such as Pseudomonas, Acinetobacter, and Vogesella exhibited higher relative abundance. As decomposition progressed into the late phase (5–14 days), the abundance of these early dominant genera significantly declined, while that of Clostridium and Paeniclostridium sharply increased, making them the dominant members of the late-stage community. This pattern was also confirmed by the stacked bar charts of genus-level composition (Figure 1B). A Sankey diagram further visualized this cross-level association, showing that in late-stage samples, the increased abundance of the phylum Firmicutes (pink flow) was almost entirely attributed to Clostridium and Paeniclostridium, whereas the early-stage Proteobacteria (blue flow) was composed of several different genera (Figure 1D). To assess the potential influence of the surrounding water on the muscle microbiome, we compared the microbial genera shared between water and muscle samples collected at day 0 and day 14. The analysis revealed distinct dynamics between the two groups (Supplementary Figure 2). In the post-mortem submersion (PS) group, the number of shared genera increased from 64 at day 0 to 111 at day 14. In contrast, the drowning (D) group exhibited a high initial number of shared genera (112) at day 0, which remained at a similarly high level (105) by day 14.

Limited impact of cause of death on microbial community diversity

3.2

To investigate whether the manner of entry into the water was a key factor influencing the microbial succession in skeletal muscle, we directly compared the community structure and diversity of the drowning (D) and post-mortem submersion (PS) groups. A principal co-ordinates analysis (PCoA) based on the Bray–Curtis distance showed that the sample points from the D and PS groups were highly intermixed in the ordination space, without forming separate, distinct clusters. The 95% confidence ellipses for the two groups also largely overlapped (Figure 2A). To statistically validate this observation, we performed a permutational multivariate analysis of variance (PERMANOVA). The results confirmed that there was no significant difference in the overall community structure (beta diversity) between the two groups (PERMANOVA, R^2^ = 0.008, p = 0.251). The R^2^ value indicated that the cause of death could only explain 0.8% of the total community variation.

Comparison of muscle microbial community diversity between the drowning (D) and post-mortem submersion (PS) groups. (A) Principal co-ordinates analysis (PCoA) plot based on Bray–Curtis dissimilarity, illustrating the overall microbial community structure. The first two principal coordinates explain 22.82 and 13.91% of the total variation, respectively. Each point represents an individual sample, colored by group (D: red; PS: blue). Ellipses indicate the 95% confidence interval for each group. (B,C) Boxplots comparing the alpha diversity between the D and PS groups at each PMSI. Community richness was estimated using the Chao1 index (B), and community diversity was estimated using the Shannon index (C). No significant differences (ns, p > 0.05) were observed between the two groups at any time point according to the Wilcoxon rank-sum test. The boxes represent the interquartile range (IQR), the central line indicates the median, and the whiskers extend to 1.5 times the IQR.

We further assessed the alpha diversity of the communities. The results showed that for both the Chao1 index, which measures species richness (Figure 2B), and the Shannon index, which measures both richness and evenness (Figure 2C), the boxplot distributions for the D and PS groups were highly similar at all examined time points (0, 1, 3, 5, 7, and 14 days). A Wilcoxon rank-sum test confirmed that at no time point was there a statistically significant difference in any alpha diversity index between the two groups (all p > 0.05). Taken together, the results from both beta and alpha diversity analyses indicate that, under the conditions of this study, the cause of death (drowning vs. post-mortem submersion) was not a major driver shaping the successional patterns of the skeletal muscle microbiome.

Stage-wise division of microbial community succession based on time

3.3

Having established that the cause of death had a limited impact, we analyzed all samples together to investigate whether the microbial community exhibited stage-wise successional characteristics. A PCoA based on the Bray–Curtis distance matrix revealed that all sample points were distributed along a clear temporal gradient, primarily driven by the first principal coordinate (PCoA1), which explained 22.8% of the community variation (Figure 3A). To further quantify this successional pattern, we divided the PMSI into an early-phase (EP; 0–3 days) and a late-phase (LP; 5–14 days). Both PCoA and non-metric multidimensional scaling (NMDS) analysis (stress = 0.1708) consistently showed that the samples from these two phases formed two statistically distinct clusters (PERMANOVA, p < 0.001), with their 95% confidence ellipses showing minimal overlap in the ordination space (Figure 3B).

Distinct microbial community structures in rat muscle at early and late post-mortem submersion intervals. (A) Principal co-ordinates analysis (PCoA) and (B) non-metric multidimensional scaling (NMDS) plots illustrating the beta-diversity of microbial communities, based on the Bray–Curtis dissimilarity metric. The first two principal coordinates in the PCoA plot explain 22.8 and 13.94% of the total variation, respectively. The NMDS plot has a stress value of 0.1708. Each point represents an individual sample. Samples are colored according to their post-mortem submersion interval (PMSI) phase: yellow for the early-phase (EP; 0–3 days) and purple for the late-phase (LP; 5–14 days). Ellipses represent the 95% confidence interval for each group. The separation between the two phases was statistically significant (PERMANOVA, p < 0.001).

PMSI prediction based on a stage-wise succession model

3.4

Based on the clear stage-wise succession of the microbial community, we evaluated the performance of several machine learning algorithms in classifying samples into the early-phase (EP, 0–3 days) and late-phase (LP, 5–14 days). We constructed four classification models: random forest (RF), logistic regression (LR), k-nearest neighbors (KNN), and support vector machine (SVM). By comparing their confusion matrices and receiver operating characteristic (ROC) curves, we found that the RF model exhibited the best classification performance (Figure 4A). Specifically, on an independent test set of 22 samples, the RF model correctly predicted 20 samples, achieving an overall accuracy of 90.9%. Its area under the ROC curve (AUC) was the highest among the four models, reaching 0.9504 (Figure 4B).

Performance evaluation of machine learning models for classifying early-phase (EP) and late-phase (LP) post-mortem submersion intervals. (A) Confusion matrices for four different classification models: random forest (RF), logistic regression (LR), k-nearest neighbors (KNN), and support vector machine (SVM). The values in each cell represent the number of samples in the test set. The color intensity corresponds to the number of samples. (B) Receiver operating characteristic (ROC) curves for the four models. The area under the curve (AUC) for each model is displayed in the legend.

After successfully establishing the stage classification, we built separate RF regression models for the EP and LP phases to accurately predict the PMSI. For the EP group (0–3 days), the model demonstrated extremely high prediction accuracy. As shown in Figure 5A, the sample points for both the training and test sets were tightly clustered around the ideal prediction line (y = x). The model’s mean absolute error (MAE) on the training set was 0.16 days, and on the independent test set, the MAE was only 0.303 days. For the more complex LP group (5–14 days), the regression model’s predictions were also robust. The scatter plot (Figure 5B) showed a strong correlation between the predicted and actual values for most samples. The MAE for this model was 0.481 days on the training set and 1.293 days on the test set.

Performance of the two-stage random forest regression models for PMSI prediction. Scatter plots showing the correlation between the observed PMSI (in days) and the PMSI predicted by the regression models for (A) the early-phase (EP; 0–3 days) and (B) the late-phase (LP; 5–14 days). Each point represents an individual sample, colored by its inclusion in either the training set (blue) or the test set (red). The solid black line represents the ideal prediction line (y = x). The mean absolute error (MAE) for both the training and test sets is displayed in the inset box for each model.

Discussion

4

This study aimed to investigate whether the microbial community of skeletal muscle, a relatively closed internal ecosystem, undergoes a predictable succession following submersion in water (PMSI) and to evaluate its utility for forensic PMSI estimation. Our central finding confirms that the rat skeletal muscle microbiome undergoes a highly predictable, stage-wise succession, enabling the construction of a two-stage machine learning model that accurately predicts PMSI up to 14 days. These results open a new avenue in forensic microbiology by highlighting the potential of internal tissues as stable sources of information for PMSI estimation.

Our research first established the overall pattern of microbial succession in skeletal muscle. Beta diversity analysis clearly showed a significant temporal gradient in community structure, which could be distinctly partitioned into an early (0–3 days) and a late (5–14 days) phase (Figure 3). This finding is consistent with previous studies on traditional sample sites like the gut and skin, which also report time-dependent microbial succession as the biological basis for the “microbial clock” (Metcalf et al., 2013; Guo et al., 2016). The novelty of our work, however, lies in systematically demonstrating this principle within the relatively “sterile” environment of skeletal muscle. This approach complements other sampling strategies aimed at estimating PMSI from skeletal remains. For instance, Rose et al. (2025) demonstrated the efficacy of using surface swabs, a valuable non-destructive method. Our focus on a well-protected internal tissue like skeletal muscle is intended to capture a more stable, endogenously-driven successional signal, thereby potentially reducing environmental stochasticity. We speculate that this phased transition marks a critical ecological turning point in the decomposition process. The early phase may be dominated by facultative anaerobes that spread via the circulatory system after death, while the late phase, characterized by increased tissue hypoxia and structural breakdown, provides an ecological niche for the proliferation of obligate anaerobes, particularly various Clostridium species.

At the taxonomic level, we observed a significant shift from Proteobacteria to Firmicutes dominance, specifically a succession from genera like Pseudomonas to Clostridium (Figure 1). This pattern not only corroborates the ecological shift mentioned above at a species level but also shares commonalities with findings from classic terrestrial decomposition models. For instance, during the bloat and active decay stages of terrestrial decomposition, obligate anaerobic Clostridium species are widely recognized as the core microbes driving tissue breakdown and gas production (Metcalf et al., 2016; Cobaugh et al., 2015). Our data suggest that even in the cold and humid conditions of an aquatic environment, the internal microenvironment of skeletal muscle (e.g., hypoxia) may still be the primary intrinsic driver of core microbial succession, with its underlying logic bearing some resemblance to that in terrestrial environments. Notably, the Sankey diagram (Figure 1D) visually demonstrates that the increase in Firmicutes in the late stage is almost entirely contributed by Clostridium and Paeniclostridium, a highly focused successional pattern that may offer promising candidates for biomarkers of long-term PMSI.

A key finding of this study is that the manner of entry into the water (drowning vs. post-mortem submersion) did not have a statistically significant effect on either the alpha or beta diversity of the skeletal muscle microbiome (Figure 2). This result may seem counterintuitive, as the drowning process theoretically introduces a large number of aquatic microbes into the respiratory and circulatory systems. A possible explanation is that skeletal muscle, as a well-vascularized but relatively closed tissue, presents a strong “ecological barrier” due to its robust anaerobic microenvironment and specific nutrient conditions. This barrier may make it difficult for typically aerobic or facultative anaerobic microbes from the water to effectively colonize and proliferate, thus failing to disrupt the overall successional trajectory dominated by endogenous bacteria. This finding has significant practical implications, suggesting that PMSI estimation models based on skeletal muscle may have greater stability than those based on surface microbes, being less affected by case-specific details such as the manner of death. This is consistent with the findings of Zhang et al. (2022) in the mouse gut, where they also observed that the cause of death had a limited impact on the succession of the gut microbiome.

The microbial exchange dynamics between water and muscle offered key insights. In the PS group, the initial sharing of 64 genera likely represents low-abundance transient and background taxa rather than significant colonization. More intriguingly, the decrease in shared taxa in the D group (from 112 to 105) suggests a process of ecological selection. We hypothesize that the selective internal microenvironment favors endogenous anaerobes that outcompete the large microbial inoculum introduced during drowning. This compelling hypothesis of competitive exclusion warrants future investigation. Crucially, the resilience of the overall successional trajectory despite this major initial perturbation strongly supports the internal muscle microbiome as a robust target for PMSI estimation.

Based on these findings, we successfully constructed and validated a high-precision, two-stage PMSI prediction model. Our findings that microbial succession on submerged bone can be used to predict PMSI are highly consistent with the foundational work in this area. Notably, Cartozzo et al. (2021a) and Cartozzo et al. (2021b) have successfully developed similar random forest models based on the microbiome of Sus scrofa bones submerged in both freshwater lake and river environments, establishing the principle that bone is a reliable substrate for long-term PMSI estimation. By first using a classification model to accurately determine the successional stage of a sample (AUC = 0.9504), and then applying a regression model optimized for that specific stage, we significantly improved the prediction accuracy (Figures 4, 5). Particularly in the early PMSI phase (0–3 days), the model’s mean absolute error was as low as 0.303 days, which is of great value for forensic applications. Even in the more complex late phase (5–14 days), the mean error of 1.293 days demonstrates its considerable potential as an auxiliary diagnostic tool. This “classify-then-regress” modeling strategy may be an effective paradigm for addressing the prediction challenges posed by the non-linear and stage-wise changes in microbial data, and is worth promoting in future forensic microbiology research (Yuan et al., 2023).

Despite the positive results of this study, several limitations exist. First, this study used an SD rat model, and there are physiological, immunological, and initial microbial compositional differences between rats and humans. Therefore, the direct applicability of the model to human cases still needs to be ultimately validated through research on other animal models (Yang et al., 2023). Second, the experiment was conducted under ambient temperature, and the lack of continuous water temperature monitoring prevented the calculation of accumulated degree-days (ADD). While the use of post-mortem submersion days served as a viable proxy in this exploratory study, we acknowledge that ADD is a more precise metric in forensic science. Future studies should prioritize the logging of water temperature to build more robust, ADD-based predictive models. Finally, this study focused primarily on the community structure of the microbiome and did not delve into the changes in its functional metabolism, the latter of which could provide deeper insights into the intrinsic mechanisms of succession and the discovery of more stable functional biomarkers.

Conclusion

5

In summary, this study is the first to systematically reveal the regular successional pattern of the rat skeletal muscle microbiome during post-mortem submersion in water, and successfully construct a high-precision two-stage PMSI prediction model. Our results not only show a novel and stable research goal for forensic microbiology but also provide a solid experimental basis and an innovative modeling approach for development of more reliable PMSI estimation tools. Future research should focus on verifying these findings in models closer to humans (e.g., pig models) and under more complex environmental conditions, and should be combined with macro-functional genomics and other methods. The ultimate goal is to transform this “microbial clock” into a powerful timing tool for forensic practice.