The Muscomorpha is one of the most species-rich brachyceran groups in Diptera, with many species serving as important disease vectors; however, its high-level phylogenetic relationships have long been controversial and unsolved. This study comparatively analyzed the characteristics of mitogenomes of 131 species that represent 18 superfamilies in Muscomorpha, in which mitogenomes of 16 species have been newly sequenced and annotated, demonstrating that their gene composition, order, AT bias, length variation, and codon usage are consistent with documented dipteran mitogenomes. The phylogenetic topologies demonstrated that the robustness of Muscomorpha and major clades within Muscomorpha are monophyletic: Cyclorrhapha, Schizophora, and Calyptratae. A clade of Empidoidea were recovered as the sister-group to Cyclorrhapha. Within Cyclorrhapha, Platypezoidea and Syrphoidea were sequentially placed as basal groups of the Cyclorrhapha. The remaining cyclorrhaph superfamilies gathered as two main clades. Ephydroidea were, in most cases, placed as the sister-group to Calyptratae. Within Calyptratae, Hippoboscoidea were sister to an assemblage of lineages composed of a Oestroid grade and Muscoidea. The Muscomorpha was proposed to originate in the early Jurassic, and the main clade diversified near the Cretaceous–Paleogeneextinction event, estimated using the MCMCtree and six fossil calibration points. The ancestral area of origin and geographic range of Muscomor- pha was deduced to be the Palaearctic region with 56.9% probability using the RASP software based on a dated tree.

Our phylogenetic analysis included published sequences from 131 muscomorpha, representing 53 families from 18 superfamilies. Multiple sequence alignment precedes matrix generation, we employed the codon-aware program MACSE v2.06 for 13 PCGs and MAFFT version 7.0 with the G-INS-i strategy for 2 rRNAs, thereafter, the 13 PCGs were subjected to trimming using Gblocks under the invertebrate mitochondrial genetic code, while the two rRNA sequences underwent trimming using trimAl v1.2rev57, subsequently, all individual alignments were concatenated into a supermatrix using the Phylosuite\_v1.2.3 platform with default settings. We constructed 4 datasets for phylogenetic analyses: (1) PCGsrRNA, the combination of 13 protein-coding genes plus two rRNA genes, resulting in a total sequence length of 12,641 nucleotides; (2) PCGs12rRNA, to mitigate substitution saturation, the third codon positions of 13 PCGs were excluded; (3) PCGs, all codon positions; and (4) AA, amino acids translated by PCGs. Before phylogenetic analyses, the substitution saturation of each codon position of the 13 mitochondrial PCGs was assessed using the index (Iss) with DAMBE v.6. the completeness of multiple sequence alignments was quantified by AliStat, and the heterogeneity of sequence was visualized using AliGROOVE v.1.08. To determine the optimal partitioning schemes and corresponding nucleotide substitution models for each dataset, we employed ModelFinder to select the best-fit substitution model for each partition in maximum likelihood (ML) analysis. To avoid the influence of heterotachous evolutionary sequences on phylogenetic inference, we used the single topology (GHOST) model in IQ-TREE, the Bayesian information criterion (BIC) and the 'greedy' algorithm were used, with branch lengths estimated as 'unlinked', to search for the best-fit scheme in the partition model. To mitigate the effects of long-branch attraction (LBA) artefacts, the posterior mean site frequency (PMSF) model was manipulated in IQ-TREE too. In the concatenated analyses, support values were assessed using the ultrafast bootstrap (UFBoot), the approximate likelihood ratio test (SH-aLRT), and a Bayes test.