Unraveling forensic timelines using molecular markers in Phormia regina maggots

Abstract

In 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 estimate the time since death. The blowfly Phormia regina is a key forensic insect, yet age estimation for older maggots in this species is particularly challenging due to the limited morphological changes in the late-stage larvae. To enhance 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. 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.

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

(A) Depiction of the sampling workflow where newly hatched maggots were transferred into a paper containing 30g chicken liver within a rearing box. The boxes were placed in a 27.5°C incubator. At each age cohort, maggots were visually scored for feeding or wandering behaviors and weighed. Created with BioRender.com. (B) 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). (C) 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.

Figure 2.. Chromosomal level genome assembly shows 6 pseudohaploid chromosomes in Phormia regina.

(A) 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 with BioRender.com. (B) 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. (C) Hi-C contact matrix heat map of the chromatin interactions and scaffolding compartments. (Blue: Chromosomes, Green: scaffolds).

Figure 3.. 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, applying a negative of log10 adjusted p value < 0.01 and log2 fold change > 1 as 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%).

Figure 4.. 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. DEGs were identified using DESeq2, applying thresholds of log₂ fold change > 1 and adjusted p-value < 0.01 (−log₁₀ scale on the y-axis). 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 (~10%).

Figure 5.. Weighted gene co-expression network analysis (WGCNA) reveals gene expression modules associated with developmental age and behavioral state 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 (blue, brown, turquoise) across samples, arranged by age. These modules exhibit coordinated gene expression patterns associated with developmental and behavioral transitions.

Figure 6.. Gene ontology enrichment analysis highlights age- and behavior-associated metabolic pathways in maggots.

This figure presents gene ontology (GO) enrichment analyses of differentially expressed genes (DEGs) across developmental timepoints (80 h, 90 h, 130 h) and behavioral states (feeding vs. wandering) in maggots. For each comparison, enriched GO terms in three categories—Biological Process, Cellular Component, and Molecular Function—were identified. Enrichment significance was determined using the Benjamini-Hochberg procedure with a false discovery rate (FDR) threshold of < 0.05. The Gene Ratio represents the proportion of DEGs associated with a specific GO term relative to the total number of input DEGs. Analyses reveal shifts in metabolic, cellular, and molecular pathways associated with age and behavior transitions. Created with BioRender.com.

Figure 7.. 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).

Figure 8.. Nine candidate genes were among all three analyses.

(A) Six candidate genes (Y5078, Y5076, Agt, Asm, Ech1 dhb4) appeared upregulated linear regression patterns. (B) Three candidate genes (Gabd, Acohc, Ivd) appeared downregulated linear regression patterns. (C) 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.