Complete Mitochondrial Genome of Piophila casei (Diptera: Piophilidae): Genome Description and Phylogenetic Implications

Piophila casei is a flesh-feeding Diptera insect that adversely affects foodstuffs, such as dry-cured ham and cheese, and decaying human and animal carcasses. However, the unknown mitochondrial genome of P. casei can provide information on its genetic structure and phylogenetic position, which is of great significance to the research on its prevention and control. Therefore, we sequenced, annotated, and analyzed the previously unknown complete mitochondrial genome of P. casei. The complete mt genome of P. casei is a typical circular DNA, 15,785 bp in length, with a high A + T content of 76.6%. It contains 13 protein-coding genes (PCG), 2 ribosomal RNA (rRNA) genes, 22 transfer RNA (tRNA) genes, and 1 control region. Phylogenetic analysis of 25 Diptera species was conducted using Bayesian and maximum likelihood methods, and their divergence times were inferred. The comparison of the mt genomes from two morphologically similar insects P. casei and Piophila megastigmata indicates a divergence time of 7.28 MYA between these species. The study provides a reference for understanding the forensic medicine, taxonomy, and genetics of P. casei.


Introduction
P. casei (Linnaeus, 1758) (Diptera: Piophilidae) belongs to the holometabolous insects, and its growth cycle includes four stages: egg, larva, pupa, and adult. The life cycle of P. casei lasts 12-30 days, and it can usually breed 7-9 generations per year [1]. P. casei exhibits a highly developed olfactory system, displaying a specific attraction to fishy and putrid odors [2,3], and has been identified as the most harmful pest to dry-cured ham. Its voracious appetite enables it to quickly consume ham, leading to the formation of cavities, darkening of the flesh, and the production of a putrid smell, ultimately resulting in significant damage to ham quality and considerable economic losses for ham manufacturers [4]. Similarly, it has been found in cheese and fish, as well as in human corpses in advanced stages of decomposition. It is widely distributed in the world and has been documented not only in Guizhou, Yunnan, and Zhejiang provinces in China but also in many parts of Europe and North America [5]. In Italy, P.
casei is used to make cheese products with a unique flavor [6]. However, ingestion of the larvae by humans can lead to intestinal myiasis, so safety remains a concern [7,8].
P. casei is generally not favored in the food industry, but it is highly valued in the fields of forensic medicine and forensic entomology. In forensic medicine, the Piophilidae family is extremely important, especially in the advanced stages of decomposition, as it can provide valuable information when determining the minimum post-mortem interval [5,9].
Piophilidae is a small family of flies that comprises several species, including P. casei and Piophila megastigmata (McAlpine 1978) [10,11]. In the early years, the high similarity in appearance between P. casei and P. megastigmata, coupled with the lack of research on P. megastigmata, made it challenging to differentiate between the two species, with P. megastigmata often being misidentified as P. casei. Until recent times, P. megastigmata was found to be more prevalent on cadavers than P. casei [12]. However, the genetic information available on P. casei is still limited, and forensic entomologists rely mostly on morphological observations rather than genetic comparisons to distinguish and compare the appearance of P. casei and P. megastigmata [13,14].
Mitochondria (mt) genes are powerful molecular markers in phylogenetic and population genetics studies. The complete mt genomes that contain more variants provide more information than partial mt gene sequences to help scientists resolve phylogenetic relationships and infer species evolution [15][16][17]. The mt genome of insects and mites generally contains 13 protein genes, 22 tRNAs, 2 rRNAs, and a control region [18][19][20]. In this study, the complete mt genome was assembled, and the basic characteristics of P. casei mt genes were studied, focusing on the comparison of the PCGs between P. casei and P. megastigmata. The known mt genome data of 25 Diptera insects from GenBank was selected for the construction of phylogenetic relations and evolution time, laying the foundation for a more in-depth study at the molecular level and providing theoretical bases and ideas for the accurate identification of forensic insects P. casei and P. megastigmata.

Sample Collection
Fifty adult specimens of P. casei were collected on 2 November 2021 from the dry-cured hams in Panxian (25 • 75.87 N; 104 • 52.57 E), Guizhou Province, China. The specimens were preserved in 100% alcohol for species identification and DNA extraction in the laboratory. Illustrations of P. casei were obtained from samples of eggs, larvae, pupae, and adults that were collected. Color images were captured using the Keyence VHX-6000 model (Keyence, Osaka, Japan) color-image analysis instrument.

DNA Extraction, Library Construction, and Sequencing
Genomic DNA was extracted from P. casei male adults using Qiagen DNeasy Blood and Tissue Kit (Qiagen, Hilden, Germany) according to the manufacturer's instructions. DNA degradation and contamination were monitored on 1% agarose gels. DNA purity was measured using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA), and DNA was quantified using Qubit DNA Assay Kit Fluorometer (Life Technologies, Carlsbad, CA, USA).
The Truseq Nano DNA HT Sample Preparation Kit (Illumina, San Diego, CA, USA) was used to generate sequencing libraries following the manufacturer's instructions, and index codes were added to assign the sequences to each sample. The DNA sample (1.5 µg DNA per sample) was fragmented by sonication to a size of 350 b. Then the DNA fragments were end-polished, A-tailed, and ligated to the full-length adapter from the sequencing kit with further PCR amplification. Last, the PCR products were purified (AM-Pure XP system) and analyzed for size distribution by an Agilent2100 Bioanalyzer. These constructed libraries were sequenced with the Illumina NovaSeq 6000 platform (BIOZE-RON, Shanghai, China), and 150 bp paired-end reads were generated with an insert size of around 350 bp [21].

Genome Assembly and Gene Annotation
In this study, a script implemented in MitoZ was utilized to filter the extracted mitochondrial gene sequences, which effectively removed reads containing numerous 'N's, low-quality reads, or PCR duplicates (i.e., pairs of identical reads) [22]. Then the mt genome was assembled in MitoZ using De novo assembly. De novo assembly is an algorithm well suited to mt mitochondrial genome assembly, which, with the assistance of other algorithms and evaluation metrics in MitoZ, enables better assembly of genomic sequences [23]. The principle is that the average sequencing depth of mitochondrial genome reads is much higher than that of the nuclear genome, and different Kmer parameters are set to achieve the best possible assembly. The Kmer parameter was finally adjusted to 33 to achieve the best assembly results. Protein-coding genes were annotated using scripts in MitoZ. tBLAST in BLAST v2.2.19 [24] was used to find candidate PCG sequences by matching sequences to the protein database sequences (MitoZ has a built-in database of insect PCG annotations); Genewise v2.2.0 [25], which produces faster and more accurate results, was used to identify each PCG for annotation. tRNAs were annotated using MiTFi [26] of the covariance model (CM) for annotation, and rRNAs were annotated using infernal-1.1.1 [27] and rRNA CMs based on an extensive manually curated alignment [28]. The complete mitogenome sequence of P. casei has been submitted to GenBank (ON204020).
To understand the taxonomic position of P. casei in Diptera, the mt genomes of 25 insects (including 2 outgroup mt genomes) were downloaded from GenBank and used to construct phylogenetic trees using both the maximum likelihood (ML) and Bayesian inference (BI) methods. The nucleotide diversity (Pi), the nonsynonymous substitution rate (Ka), and the synonymous substitution rate (Ks) were calculated using Launch DnaSP6 v6.12.3 [30], and the results were plotted using GraphPad Prism v8.0.2. Base sequence substitution saturation was analyzed and plotted using DAMBE v7.3.11 [31].
The 22 tRNA genes, 13 PCGs, and 2 rRNA genes of each of the 25 insects were extracted from NCBI and imported into PhyloSuite v1.2.2 [32]. The sequences were trimmed using Gblocks v0.91 b [33] to remove redundant codons and then concatenated using concatenate sequence to form a single concatenated sequence. All sequences were then aligned using MAFFT v7.313 and optimized using MACSE v2.0.1 [34]. The secondary structure of 22 tRNAs was predicted with ARWEN and compared manually. The MFE (minimum free energy) structures of two ribosomal RNA genes (rrnL and rrnS) were obtained by RNAfold [35] prediction. For the tandem sequences, the optimal GTR + F + I + G4 and GTR + F + G4 molecular phylogenetic model was found using ModelFinder [36] and imported into MrBayes v3.2.6 [37] and IQ-TREE v1.6.8 [38] for BI and ML phylogenetic tree construction with 104 and 2 × 106 bootstraps, respectively. iTol (https://itol.embl.de/itol. cgi (accessed on 10 September 2022)) was used to visualize the phylogenetic trees, and then the developmental trees constructed by the two methods were compared.

Divergence Time Estimate
The divergence time was estimated from the BI phylogenetic trees with timetree5 (http: //www.timetree.org/ (accessed on 10 September 2022)) and BEAST v1.10.4 [39] according to the related literature and using age 70 MYA as the divergence time calibration [40]. The estimation of divergence models for subfamilies was based on complete mt genomes according to the strict clock log-normal model in BEAST. A calibrated Yule model was used and spliced using the GTR + F + I + G4 and GTR + F + G4 models. After confirming the convergence of the chains with Trancer v1.7.2 [41], the first 10% of generations were burn-in as ageing every 1000 generations sampled at 10 7 conditions. We summarized the subsample trees in a maximum clade credibility tree with mean heights using Tree Annotator v1.10.4. The mean heights and 95% highest probability density (95% HPD) were displayed in Figtree v1.4.3. The labeling and photo-compositing of the images in this article were made using Adobe Photoshop 2022 and Adobe Illustrator 2022.

General Features of P. casei mt Genome
The complete mitochondrial genome of P. casei is 15,785 bp long (GenBank No. ON204020) and contains 13 PCGs (Protein-coding genes) (11,206 bp), 22 tRNA genes (1467 bp), 2 rRNA genes (2107 bp), and 1 control region (744 bp). The circular map of the P. casei mt genome is shown in Figure 1a and is generally consistent with the mt genomes of eight other Diptera, such as P. megastigmata (15, (Figure 1f). The difference in mt genome length between these five species is mainly due to the differences in the size of the control region (from a minimum of 388bp to a maximum of 1354 bp), as previously observed between Diptera species [42]. In Diptera, the mt genomes of some species, such as Zeugodacus caudatus (15,311 bp) (Figure 1d  The mitochondrial genomes of P. casei and P. megastigmata exhibit a high degree of conservation, as evidenced by the presence of 37 genes in each genome, as well as the absence of any inversion events or swapping of gene coding directions. In terms of ATCG content and AT-skew, CG-skew, the mt genome of P. megastigmata is only slightly different than that of P. casei ( Table 1). The difference between the two is more pronounced in the AT% and AT-skew of the tRNA gene, indicating differences in amino acid content between the two species. The complete P. casei mt genome consisted of 39.5% A, 37.1% T, 13.8% C, and 9.6% G. The total A + T content was 76.6%, which is relatively high among Tephritidae species. The mt genome had a positive AT skew (0.03) and a negative CG skew (−0.18); according to the previous reports, a bias against the use of Gs in all strands is characteristic of the metazoan mitochondrial genome [43]. Four PCGs (nad5, nad4, nad4l, nad1), eight tRNA genes (trnQ, trnC, trnY, trnF, trnH, trnP, trnL, trnV), and two rRNA genes (rrn12, rrn16) were found on the minority strand (N strand), while the remaining genes were on the majority strand (J strand) ( Table 2). Ten genes overlap on the mt genome of P. casei, with a total of thirty-four overlapping bases. There are 14 gene spacers, with the largest spacer (16 bp) occurring between rrn12 and the control region.
In order to compare the AT% variations among different species within the same family, nine species in Tephritoidea were selected [41]. Overall, the AT content of each of their genes was found to be similar, with significant differences only in the proportion of codons used in the first, second, or third codons. The first codon position in both P. casei and P. megastigmata was found to have higher AT% compared with several other species in the Tephritoidea (Figure 2).  In terms of ATCG content and AT-skew, CG-skew, the mt genome of P. megastigmata is only slightly different than that of P. casei ( Table 1). The difference between the two is more pronounced in the AT% and AT-skew of the tRNA gene, indicating differences in amino acid content between the two species.
The complete P. casei mt genome consisted of 39.5% A, 37.1% T, 13.8% C, and 9.6% G. The total A + T content was 76.6%, which is relatively high among Tephritidae species. The mt genome had a positive AT skew (0.03) and a negative CG skew (−0.18); according to the previous reports, a bias against the use of Gs in all strands is characteristic of the metazoan mitochondrial genome [43]. Four PCGs (nad5, nad4, nad4l, nad1), eight tRNA genes (trnQ, trnC, trnY, trnF, trnH, trnP, trnL, trnV), and two rRNA genes (rrn12, rrn16) were found on the minority strand (N strand), while the remaining genes were on the majority strand (J strand) ( Table 2). Ten genes overlap on the mt genome of P. casei, with a total of thirty-four overlapping bases. There are 14 gene spacers, with the largest spacer (16 bp) occurring between rrn12 and the control region.
In order to compare the AT% variations among different species within the same family, nine species in Tephritoidea were selected [41]. Overall, the AT content of each of their genes was found to be similar, with significant differences only in the proportion of codons used in the first, second, or third codons. The first codon position in both P. casei and P. megastigmata was found to have higher AT% compared with several other species in the Tephritoidea (Figure 2).
The nucleotide diversity (Pi) was compared for the PCGs of 26 species (Figure 3a), and the strength of polymorphisms was used to characterize the PCGs of 26 species. Overall, the Pi values ranged from 0.153 to 0.269. The nad2 (Pi = 0.269) gene had the greatest variability among these genes, followed by nad6 (Pi = 0.259), atp8 (Pi = 0.255), nad4l (Pi = 0.202), and, finally, cox1 (Pi = 0.153) with the least variability. It shows that cox1 has high genetic stability among dipterans, while nad2 has weak stability and is suitable for the study of the evolution of variation [45]. The Pi values of P. casei and P. megastigmata were also compared individually in each species (Figure 3b): cox1 (Pi = 0.018) demonstrated the least variability, nad3 (Pi = 0.120) showed the greatest variability, and atp8 (Pi = 0.025) had a low degree of variability, unlike the results of the multispecies comparison.  The frequency of non-synonymous mutations (Ka), the frequency of synonymous mutations (Ks), and the ratio of the two (Ka/Ks) were calculated for 13 PCGs of 26 species and used to understand the rate at which synonymous and non-synonymous mutations occur (Figure 3d). The Ka/Ks values determine whether there is selection pressure on the PCGs and also the degree of conservation of the coding gene. The closer the Ka/Ks is to 1, the lower the selection pressure is on that gene. The Ka/Ks values range from 0.0561 to 0.5119, and the order of Ka/Ks values for PCGs from the smallest to the largest is as follows: cox1, cox3, cytb, cox2, atp6, nad1, nad3, nad5, nad4, nad2, nad6, nad4l, atp8. Among them, the genes of the cytochrome oxidase (cox) family cox1, cox2, cox3, and cytb have smaller Ka/Ks values, indicating that they face a stronger purifying selection, a slower evolutionary rate, especially cox1, which, because of the lowest Ka/Ks value, has been The frequency of non-synonymous mutations (Ka), the frequency of synonymous mutations (Ks), and the ratio of the two (Ka/Ks) were calculated for 13 PCGs of 26 species and used to understand the rate at which synonymous and non-synonymous mutations occur (Figure 3d). The Ka/Ks values determine whether there is selection pressure on the PCGs and also the degree of conservation of the coding gene. The closer the Ka/Ks is to 1, the lower the selection pressure is on that gene. The Ka/Ks values range from 0.0561 to 0.5119, and the order of Ka/Ks values for PCGs from the smallest to the largest is as follows: cox1, cox3, cytb, cox2, atp6, nad1, nad3, nad5, nad4, nad2, nad6, nad4l, atp8. Among them, the genes of the cytochrome oxidase (cox) family cox1, cox2, cox3, and cytb have smaller Ka/Ks values, indicating that they face a stronger purifying selection, a slower evolutionary rate, especially cox1, which, because of the lowest Ka/Ks value, has been subjected to the highest purity selection. This may be why it has been used for the study of evolutionary biology as the molecular marker. Combined with Pi and Ka/Ks values, the evolutionary rate and variability of cox1 are very low. Based on the computational information, the molecular weight (Da), isoelectric points, instability index, aliphatic index, GRAVY, and other indices of the proteins were predicted ( Table 2).
The Pi and Ka/Ks values of atp8 are high compared with other genes. In contrast, another atp synthase atp6 shows a much slower evolutionary rate. The genes atp6 and atp8 are arranged in close proximity and overlap by 7 bp, indicating that ap6 and atp8 did not co-evolve as a gene cluster during the evolutionary process. The Pi of atp8 within the dipteran range showed a large difference in Pi values between the two species, and the difference between the two results is a better indication that the two species share similarities in an evolutionary direction.
The 13 PCGs were further analyzed by calculating the effective codon number (ENc) (Figure 4b), with values closer to 20 indicating a higher expression and a greater codon preference. The GC of the silent third codon posit (GC3s) value was calculated to reflect the probability of synonymous codon use. Plotting the ENc-plot according to the formula Enc = 2 + GC3 + 29/(GC32 + (1 − GC3)2) with the ENc value as the vertical axis and GC3 as the horizontal axis revealed that most of the codons were distributed around the standard curve (Figure 4b). The results suggest that the formation of codon bias in the P. casei mt genome may be related to mutations alone and is less affected by other factors (e.g., natural selection) [46]. In comparison, P. megastigmata is more subject to external factors that produce codon bias (Figure 4f).  The number of synonymous codons (L_sym) and the total number of amino acids (L_aa) were analyzed [47]. Based on the computational information, the molecular weight (Da), isoelectric points, instability index, aliphatic index, GRAVY, and other indices of the proteins were predicted. The results show that cox1, cox2, cox3, and cytb are structurally stable and can be used for species classification and mutation-related studies (Table 3), which is consistent with the Pi values as well as with the Ka/Ks values analysis [48].  The number of synonymous codons (L_sym) and the total number of amino acids (L_aa) were analyzed [47]. Based on the computational information, the molecular weight (Da), isoelectric points, instability index, aliphatic index, GRAVY, and other indices of the proteins were predicted. The results show that cox1, cox2, cox3, and cytb are structurally stable and can be used for species classification and mutation-related studies (Table 3), which is consistent with the Pi values as well as with the Ka/Ks values analysis [48]. In P. casei, the five codons with the highest relative synonymous codon usage rate (RSCU) [49] are UUA (4.79), CGA (2.97), UCU (2.74), GCU (2.35), and UCA (2.31) (Figure 4a,d). The RSCU values of these codons are all greater than one, indicating that they are more frequently used in encoding amino acids. The RSCU of P. casei and P. megastigmata show a similar bias, with only minor differences in the use of a few codons, such as CGC, CUC, and ACG. As with other Diptera, the mt genome of Piophilidae is more biased toward the use of amino acids encoded by codons with A or T in the third position [42]. The amino acid use frequencies of PCGs in the mt genomes of P. casei and P. megastigmata were counted and showed that the number of 20 amino acids varied, but both had the highest frequency of leucine (Leu) usage, followed by isoleucine (Ile), serine (Ser), and glycine (Gly) (Figure 4c,g). Overall, there is still a predominance of nonpolar amino acids (hydrophobic) and very little amount of basic amino acids (Argue, Lysol, His), with no bias toward acidic (Asp, Glu). By comparing the amino acid content with the codon, it was found that the codon content with the third codon being A or G had an effect on the amino acid content. The small difference in the frequency of amino acid use between the two species suggests that there is also a bias toward the use of amino acids.

Transfer and Ribosomal RNA Genes
The complete mt genome of P. casei included 22 tRNA genes and 2 rRNA genes. The 22 tRNAs had a length between 65 bp and 72 bp. A total of 21 typical trilobal secondary structures were predicted. The tRNA gene trnS was the only one missing the D-arm ( Figure 5). It is very common in the fly family, and the absence of this arm does not affect the function of the tRNA genes such as Chrysomya chloropyga [43]. The predicted secondary structure maps of rrnL and rrnS by energy-minimization (MFE) are in general agreement with the secondary structures proposed for other insects and provide a reference for subsequent studies ( Figure 6) [19,50,51]. The features of tRNAs were conserved between the P. casei and P. megastigmata genomes. The 16s-rRNA gene is 1320 bp and is located between trnL and trnV. The 12s-rRNA gene is 787 bp and is located between trnV and the control region. The complete mt genome of P. casei has a 744 bp control region, which is responsible for regulating DNA replication and transcription. The A + T content of this region is 92.61% between rrnS and trnL.

Phylogenetic Analysis
We used the complete mt genome sequence of P. casei for phylogenetic analysis, starting with a tree-building assessment of all species' gene sequences (PCGs and RNAs). The saturation plots in Figure 7 show the relationships between general time reversible (GTR) distance and transition/transversion rates. No significant substitution saturation was detected for any codon positions of the 13 PCGs and RNAs in either the symmetrical or asymmetrical topology tests (Iss (0.2701, 0.3261) < Iss.c (0.8452, 0.8046), p < 0.001) [52]. The Iss values were all lower than the Iss.c values, indicating that the results of the tree construction using the complete mt genome of 26 species were meaningful for phylogenetic tree construction. Therefore, we subsequently used the complete mt genomes of 26 species for phylogenetic tree construction and analysis.

Phylogenetic Analysis
We used the complete mt genome sequence of P. casei for phylogenetic analysis, starting with a tree-building assessment of all species' gene sequences (PCGs and RNAs). The saturation plots in Figure 7 show the relationships between general time reversible (GTR) distance and transition/transversion rates. No significant substitution saturation was detected for any codon positions of the 13 PCGs and RNAs in either the symmetrical or asymmetrical topology tests (Iss (0.2701, 0.3261) < Iss.c (0.8452, 0.8046), p < 0.001) [52].
The Iss values were all lower than the Iss.c values, indicating that the results of the tree construction using the complete mt genome of 26 species were meaningful for phylogenetic tree construction. Therefore, we subsequently used the complete mt genomes of 26 species for phylogenetic tree construction and analysis. In this study, phylogenetic trees were constructed using 12 PCG genes, 2 rRNA genes, and 22 tRNA genes from 26 species using different best-fit substitution models, the ML algorithm (Figure 8), and the BI algorithm ( Figure 9). Two common outgroups (Culicoidea, A. oryzalimnetes, NC_030715, and Chironomoidea, Simulium variegatum, NC_033348) were used. Both ML and BI phylogenetic analyses produced a similar topology with slightly different node support values. The node support was generally higher for the BI tree than for the ML tree, which is common in studies of taxa [53][54][55]. High support was obtained throughout the phylogenetic tree, with Tephritoidea and Sciomyzoidea appearing distinctly separated from other families, and the abdomen of species within Tephritoidea is consistent with Hoi-Sen's study [56].  In this study, phylogenetic trees were constructed using 12 PCG genes, 2 rRNA genes, and 22 tRNA genes from 26 species using different best-fit substitution models, the ML algorithm (Figure 8), and the BI algorithm ( Figure 9). Two common outgroups (Culicoidea, A. oryzalimnetes, NC_030715, and Chironomoidea, Simulium variegatum, NC_033348) were used. Both ML and BI phylogenetic analyses produced a similar topology with slightly different node support values. The node support was generally higher for the BI tree than for the ML tree, which is common in studies of taxa [53][54][55]. High support was obtained throughout the phylogenetic tree, with Tephritoidea and Sciomyzoidea appearing distinctly separated from other families, and the abdomen of species within Tephritoidea is consistent with Hoi-Sen's study [56].

Figure7. Substitution saturation plots of mt genomic PCGs and RNAs for 26 species.
Plots in blue and green indicate transition and transversion, respectively.
In this study, phylogenetic trees were constructed using 12 PCG genes, 2 rRNA genes, and 22 tRNA genes from 26 species using different best-fit substitution models, the ML algorithm (Figure 8), and the BI algorithm ( Figure 9). Two common outgroups (Culicoidea, A. oryzalimnetes, NC_030715, and Chironomoidea, Simulium variegatum, NC_033348) were used. Both ML and BI phylogenetic analyses produced a similar topology with slightly different node support values. The node support was generally higher for the BI tree than for the ML tree, which is common in studies of taxa [53][54][55]. High support was obtained throughout the phylogenetic tree, with Tephritoidea and Sciomyzoidea appearing distinctly separated from other families, and the abdomen of species within Tephritoidea is consistent with Hoi-Sen's study [56].   We did not use the older calibration points (238.5-295.4 MYA) and instead obtained divergence time data that is more realistic [40]. Combined with the timing of divergence, the evolutionary history of Oestroidea species and Piophilidae species differ significantly, although they share similar feeding habits, with their ancestors diverging as early as 69.85 MYA (95% HPD: 68.69-70.89 MYA) around the late Cretaceous periods. The Piophilidae of the Tephritoidea diverged from other plant-feeding species at 46.51 MYA (95% HPD: 37.16-55.81 MYA) to produce two different diets, despite sharing similar morphological characters.
The results of the phylogenetic tree show that P. casei is very closely related to P. megastigmata; the mt genetic distance within the same family of 26 Piophilidae species is 0.062. Judging by the evolutionary branch lengths, the P. megastigmata has longer branch lengths, reflecting its greater degree of genetic variation and a higher degree of evolution. The divergence between P. casei and P. megastigmata was inferred at 7.28 MYA (95% HPD: 3.64-18.20 MYA) around the late Cenozoic Tertiary during the Miocene [42]. These results may provide a valuable basis for the study of the evolution of the P. casei gene family, population transmission, and biodefense potential [57][58][59].

Conclusions
Previous studies have attempted to describe the phylogenetic relationships of P. casei and to identify it taxonomically based on morphological characters or nucleotide fragments. In this study, we report the first complete mt genome of P. casei and document the phylogenetic relationships of P. casei with the mt genomes of other 26 fly species. This research provides the evidence for the accurate identification of two morphologically closely related forensic insects P. casei and P. megastigmata, and a theoretical basis for a more in-depth study of their genes. Unfortunately, the sample size of the species used in the article is limited, and further investigations are needed to compare the mitochondrial genome among families of Diptera as a whole, as well as to explore the species divergence time. A deeper exploration of the P. casei genome also merits further investigation. We did not use the older calibration points (238.5-295.4 MYA) and instead obtained divergence time data that is more realistic [40]. Combined with the timing of divergence, the evolutionary history of Oestroidea species and Piophilidae species differ significantly, although they share similar feeding habits, with their ancestors diverging as early as 69.85 MYA (95% HPD: 68.69-70.89 MYA) around the late Cretaceous periods. The Piophilidae of the Tephritoidea diverged from other plant-feeding species at 46.51 MYA (95% HPD: 37.16-55.81 MYA) to produce two different diets, despite sharing similar morphological characters.
The results of the phylogenetic tree show that P. casei is very closely related to P. megastigmata; the mt genetic distance within the same family of 26 Piophilidae species is 0.062. Judging by the evolutionary branch lengths, the P. megastigmata has longer branch lengths, reflecting its greater degree of genetic variation and a higher degree of evolution. The divergence between P. casei and P. megastigmata was inferred at 7.28 MYA (95% HPD: 3.64-18.20 MYA) around the late Cenozoic Tertiary during the Miocene [42]. These results may provide a valuable basis for the study of the evolution of the P. casei gene family, population transmission, and biodefense potential [57][58][59].

Conclusions
Previous studies have attempted to describe the phylogenetic relationships of P. casei and to identify it taxonomically based on morphological characters or nucleotide fragments. In this study, we report the first complete mt genome of P. casei and document the phylogenetic relationships of P. casei with the mt genomes of other 26 fly species. This research provides the evidence for the accurate identification of two morphologically closely related forensic insects P. casei and P. megastigmata, and a theoretical basis for a more in-depth study of their genes. Unfortunately, the sample size of the species used in the article is limited, and further investigations are needed to compare the mitochondrial genome among families of Diptera as a whole, as well as to explore the species divergence time. A deeper exploration of the P. casei genome also merits further investigation.