Unraveling forensic timelines using molecular markers in Phormia regina maggots

Abstract

In the medico-legal application of forensic entomology, estimating the time of death is critical and traditionally relies on changes in observable traits of carrion feeding insect larvae. Traits such as size, weight, and morphology can be used to predict the insect specimen age and help define the minimum time since death. The blowfly Phormia regina Meigen (Diptera: Calliphoridae) is a key forensic insect, yet age estimation for older maggots in this and other carrion-feeding species is particularly challenging due to the limited morphological changes in the late-stage larvae. To enhance age-estimation precision, we employed transcriptomic profiling on blowfly maggots, aiming to identify genes as markers for time of death estimation. Our study characterized maggot development, reinforcing that weight and behavior cannot precisely determine age between 100 and 130 hours at 27.5 °C. We built a chromosomal scale annotated genome, establishing a reliable database for uncovering transcriptomic signatures during larval development. Applying differential gene expression analyses, weighted gene co-expression network analysis, and the generalized linear model, we identified nine candidate genes (y5078, y5076, agt2, ech1, dhb4, asm, gabd, acohc, ivd) that delineate the age of otherwise indeterminate maggots. This research introduces a molecular approach to address a longstanding problem in forensic entomology and promises to increase precision in determining the time of death at a crime scene.

Fig 1. Tracking maggot weight and behavior from 70 hours old to 130 hours old.

(A) Scatter plot of individually weighed maggots scored for feeding and wandering from 70 to 130 h. The straight line is the mean weight at each age cohort. The graph shows change in larval weight over four replicate experiments (n = 823). Brown-Forsythe and Welch ANOVA tests were performed on the weight of aging cohorts (P < 0.0001). (B) The mean weight of feeding (yellow) and wandering (blue) maggots from 70 to 130 h. Population transition between feeding and wandering maggots occurred between 100 ~ 120 h. The graphs showed the number of maggots scored as feeding or wandering (a) over four replicate experiments (n = 823). The error bars indicate the standard error of each cohort.

Fig 2. Genome assembly pipeline workflow.

PacBio HiFi reads were assembled and duplicate sequences were purged. After the filtration of bacterial contaminant contigs from the assembly, the Hi-C sequences were mapped to the clean contigs for generating contact matrix and contact heat map. Created in BioRender. Lin, S. (2025) [32].

Fig 3. Chromosomal level genome assembly shows 6 pseudohaploid chromosomes in Phormia regina.

(A) Snail plots summarized the BUSCO, AT & GC content and scaffold statistics of P. regina genome assembly without proteobacteria contigs. Red segments represented the size of scaffold; dark and light orange segments indicated N50 and N90 values; central light gray indicated the cumulative scaffold counts under order of magnitude; blue and light blue indicate the proportion of AT & GC content. Upper right corner showed BUSCO completeness score against the Diptera database. (B) Hi-C contact matrix heat map of the chromatin interactions and scaffolding compartments. These six pseudohaploid chromosomes were clearly delineated by blue boundaries, with red signals indicating high-frequency chromatin interactions within chromosomes. (Blue: Chromosomes, Green: scaffolds).

Fig 4. Age-dependent transcriptional changes associated with the transition from feeding to wandering in P. regina maggots.

(A–F) Volcano plots showing differentially expressed genes (DEGs) between feeding and wandering larvae across different age cohorts: (A) 90 h, (B) 100 h, (C) 110 h, (D) 120 h, (E) 130 h, and (F) combined age groups. Differential expression was determined using DESeq2, with thresholds set at fold change > 2 (|log₂ fold change| > 1) and adjusted p-value < 0.01 (−log₁₀ padj > 2). In the volcano plot, the x-axis represents log₂ fold change, and the y-axis represents −log₁₀ adjusted p-value, with dashed lines indicating these thresholds. Genes are color-coded: upregulated (magenta), downregulated (cyan), and nonsignificant (gray). Notable genes such as 4ebp, agt2, and tsal are labeled in panels D and F. Statistical significance was assessed using the Wald test. NS denotes genes without significant differences in expression (~10%).

Fig 5. cstB is differentially expressed between the feeding and wandering stage maggots in both Chrysomya rufifacies and Phormia regina.

(A) Principal component analysis (PCA) of gene expression profiles between feeding and wandering third instar maggots. PCA was performed using normalized gene expression data across all treatments to assess transcriptional differences between feeding (orange circles) and wandering (blue triangles) stages. PC1 and PC2 accounted for 91.19% and 8.81% of the total variance, respectively. (B) Venn diagram comparing the gene expression changes associated with feeding and wandering in third instar maggots from Lucilia sericata [24,25] and Chrysomya rufifacies [26], and Phormia regina. The numbers indicate the total number of differentially expressed genes reported in each study.

Fig 6. Age-specific gene expression signatures in P. regina maggots across developmental cohorts.

(A–F) Volcano plots showing differentially expressed genes (DEGs) in each larval age group compared against all other age cohorts: (A) 80 h, (B) 90 h, (C) 100 h, (D) 110 h, (E) 120 h, and (F) 130 h. Differential expression was determined using DESeq2, with thresholds set at fold change > 2 (|log₂ fold change| > 1) and adjusted p-value < 0.01 (−log₁₀ padj > 2). In the volcano plot, the x-axis represents log₂ fold change and the y-axis represents −log₁₀ adjusted p-value, with dashed lines indicating these thresholds. Genes are color-coded as follows: upregulated (magenta), downregulated (cyan), and nonsignificant (gray). Notable differentially expressed genes such as pero, asm, armet, and hyou1 are labeled in panel B. Statistical significance was assessed using the Wald test. NS denotes genes without significant differences in expression.

Fig 7. Weighted gene co-expression network analysis (WGCNA) reveals gene expression modules associated with age in maggots.

(A,) Heatmaps of module–trait correlations showing representation of normalized reads of sample clusters (column dendrograms) by developmental age. Color scale indicates the strength of correlation between module eigengenes. (B) Expression profiles of genes within selected WGCNA modules (magenta, orange, black, cyan) across samples, arranged by age. These modules exhibit coordinated gene expression patterns associated with maggot age.

Fig 8. Age-associated candidate P. regina transcripts identified by linear regression analysis.

Linear regression was used to identify transcripts with consistent changes in expression across developmental time points (80 h to 130 h). A total of 59 transcripts were classified as upregulated (A), showing increased expression over time with positive regression slopes (β₁ > 1) and strong model fit (R² ≥ 0.91). On the other hand, 45 transcripts were classified as downregulated (B), exhibiting decreased expression with negative slopes (β₁ < –1) and R² ≥ 0.91. Gene lists are ranked by correlation coefficient, and expression is reported in transcripts per million (TPM).

Fig 9. Nine genes were candidates for estimating maggot age across all three statistical analyses.

(A) This Venn diagram displays a total of nine candidate genes that were consistently identified in all three analyses. GLM (General Linear Regression), WGCNA (Weighted Gene Co-Expression Network Analysis), and DE (Differentially Expressed Genes) analysis. The diagram illustrates the numbers of expressed transcripts associated with these three analytical approaches. Within those nine candidate genes, (B) six candidate genes (y5078, y5076, agt, asm, ech1 dhb4) appeared upregulated linear regression patterns. (C) Three candidate genes (gabd, acohc, ivd) appeared downregulated linear regression patterns.