All anthomyiids were captured by malaise traps in the Baihua Mountain (39°50′11.04″N, 115°34′41.52″E) and Dalaoling National Natural Reserve (31°4′35.6″N and 110°56′11.6″E), from 2017 to 2019 in China. All experimental materials were preserved in absolute ethanol and cryopreserved at –20°C until further processing in the Museum of Beijing Forestry University (BFU), Beijing, China. Specimens of Anthomyiidae were initially identified by Mingfu Wang using available taxonomic keys (Xue and Chao 1998), and identifications were confirmed using DNA barcodes (cox1) obtained from the assembled mitogenomes held in public databases (i.e., BOLD, NCBI) and confirmed by BLAST search to the genus level (Michelsen 2011). All Anthomyiidae mitogenome data from NCBI were downloaded and employed in comparative mitogenomic analyses with the 18 new mitogenomes in this study (Table 1).

We used the DNeasy Blood and Tissue kit (Qiagen, Hilden, Germany) according to manufacturer’s protocol for DNA extracted from individual adult flies. Qubit 3.0 was used to quantify the concentration of the DNA samples. To enhance sequencing efficiency and minimize resource waste, hybrid libraries were adopted (Gillett et al. 2014). Subsequently, the genomic DNA was pooled and then sequenced on the Illumina Novaseq 6000 platform (PE150, Illumina, San Diego, CA). Raw reads were trimmed using Trimmomatic (Bolger et al. 2014), with each library yielding approximately 5 Gb of clean data. These were assembled de novo using IDBA-1.1.1 (Peng et al. 2012). To identify mitogenomes, two sequence fragments of mtDNA (cox1 and cytb) (Crampton-Platt et al. 2015; Yan et al. 2019) were amplified as bait sequences to acquire the best-fitting mitochondrial scaffolds using Basic Local Alignment Search Tool (BLAST) with a similarity threshold of 98% (Altschul et al. 1990). The 13 protein-coding genes (PCGs) and two ribosomal RNA genes (rRNAs) were annotated using Geneious v2020.0.2 by alignment to other reported Calyptratae flies for each orthologous gene (Kearse et al. 2012). Positional annotation of 22 transfer RNA genes (tRNAs) was achieved using the online MITOS tool (Bernt et al. 2013). Complete mitochondrial genomes were submitted to NCBI under the accession numbers of OP616784-OP616801. The associated SRA, BioProject, and Bio-Sample numbers are SRR25463435-SRR25463439, PRJNA1000204, and SAMN36763070-SAMN36763087, respectively.

Sequence comparisons were carried out in PhyloSuite software (Zhang et al. 2020) to estimate nucleotide composition and relative synonymous codon usage (RSCU) among the 18 newly sequenced mitochondrial genomes. Base composition skewness analysis was calculated on all available anthomyiid mitogenomes using the specific formulas: AT-skew = (A - T) / (A + T) and GC-skew = (G - C) / (G + C) (Perna et al. 1995). Nucleotide divergence (Pi) value of three subfamilies was computed through DnaSP v6. (Rozas et al. 2017). Additionally, the ratios of Ka (nonsynonymous substitutions)/Ks (synonymous substitutions) based on 13 aligned PCGs were also measured with DnaSP v6 to compare substitution rate (Rozas et al. 2017). The Kimura 2-parameter model in MEGA 5 was used in calculations of mean genetic distances among the three subfamilies (Tamura et al. 2011). DAMBE 7.0 was applied to assess the substitution saturation (Iss) of each codon position based on all PCGs under the GTR model (Xia 2018).

The 29 complete mitogenomes from three subfamilies of Anthomyiidae were chosen to construct the phylogenetic tree, including 18 new mitogenomes documented in this study. Eight outgroups were selected to represent seven outgroup families of Diptera, with the placed between Drosophilidae (Drosophila mercatorum) and all other sampled flies. Phylogenetic relationships were inferred from analyses of a dataset of the 13 mitochondrial PCGs. To construct this dataset, each PCG of 37 mitogenomes was individually aligned using MAFFT (Katoh and Standley 2013). The optimal partitioning schemes and best-fitting model for each PCG were obtained by PartitionFinder 2 (Lanfear et al. 2017). Phylogenetic analyses (ML and BI) were performed on a concatenated 13 PCG dataset using the online CIPRES Science Gateway (Miller et al. 2010). For ML analysis, the node support values were inferred by ultrafast bootstrap resampling (BP) with 1000 replicates in IQ-TREE. Two separate Markov chain Monte Carlo (MCMC) chains were carried out for BI analyses, spanning 10 million generations simultaneously, with sampling occurring every 1000 iterations. In Bayesian analyses, posterior probabilities (PPs) were calculated after discarding the initial 25% samples as burn-in. Convergence was assessed by confirming that the average standard deviation of split frequencies was less than 0.01 in MrBayes 3.2.6 and effective sample size (ESS) was greater than 200 in Tracer (Ronquist et al. 2012; Rambaut et al. 2018). Phylograms were modified and visualized using FigTree v 1.4.

Our newly sequenced mitogenomes of Anthomyiidae show some variation in genome size, ranging from 15,635 bp to 21,098 bp in length. They are compact circular, double-stranded molecules, and are composed of the core 37 genes and a control region. The majority strand (J-strand) encoded 23 genes (9 PCGs, and 14 tRNAs), while the remaining genes were transcribed on the minority strand (N-strand) (4 PCGs, 8 tRNAs, and 2 rRNAs). All newly sequenced Anthomyiidae mitogenomes were conserved in gene order and orientation, consistent with previously published Muscoidea mitogenomes (Ren et al. 2019; Oliveira et al. 2008; Li et al. 2016). All PCGs began with a typical start codon (ATN), except for the cox1 initiated with TCG. Furthermore, most PCGs ended with the termination codons TAA/TAG, the occurrence of the TAA is more frequently observed than TAG, while three PCGs (cox2, nad4 and nad5) terminated with T, which is a common phenomenon in Calyptratae (Ren et al. 2019; Oliveira et al. 2008; Li et al. 2016; Li et al. 2020; Yan et al. 2021b; Zhao et al. 2013).

Relative synonymous codon usage (RSCU) values of all three subfamilies are illustrated in Figure 1. All PCGs encoded 22 standard amino acids. The most continually encoded amino acids in Anthomyiidae are Ala, Arg, Gly, Leu2, Pro, Ser2, Thr and Val, with Ser2 in possession of the highest value across all three subfamilies. UUA (Leu2) is the most frequently utilized codon among all sampled anthomyiids.

The AT content of PCGs from all anthomyiids varies between 76.7% (Anthomyiinae) and 77.7% (Pegomyinae) with a mean value of 77%, whereas the basic composition of all PCGs is homogeneous. Third codon positions possess a substantially higher AT content than first and second codon positions, according to analyses of the average base composition at each codon position (Table 2). By contrasting the AT content of each PCG across all Anthomyiidae, nad6 (84.16%) shows the highest mean value, followed by atp8 (83.24%) and nad4L (82.16%). On the other hand, the average value of cox1 (70.85%) and cox3 (71.76%) is the lowest. Across all the 13 PCGs of Anthomyiidae, the AT-skew is negative, but is highest in Anthomyiinae (-0.152) and lowest in Myopininae (-0.157); whereas the GC-skew is positive, ranging from 0.028 to 0.038, and remains consistent across both Anthomyiinae and Myopininae (Fig. 2). This pattern resembles that observed in most Calyptratae mitogenomes, and suggests that strand bias may be a dissymmetric mutation occurring during DNA replication (Ren et al. 2019; Oliveira et al. 2008; Li et al. 2016; Li et al. 2020; Yan et al. 2021b; Zhao et al. 2013).

Additionally, saturation plots reveal that only the third codon locations in all of the PCGs exhibit notable heterogeneity, implying that levels of heterogeneity in the PCG123 are exceedingly low (Fig. 3). There is also no evidence of saturation from sequence comparisons including the more freely evolving third codon position (Yan et al. 2019), thereby supporting the suitability of nucleotide analyses for phylogenetic reconstruction at this level.

To explore sequence evolution among the 13 PCGs sampled in anthomyiids, the values of Ka, Ks, and Ka/Ks (ω) for each PCG were computed, respectively (Fig. 4). The range of Ka for the typical gene-specific substitution rates was 0.034 (cox2) to 0.375 (cox1). The gene (nad1) showed the highest evolutionary rate (ω = 1.001) of all the PCGs, indicating that it is likely undergoing positive or relaxed selection pressure. In contrast, cox2 displayed the lowest value (ω = 0.108), reflecting that may be subject to strong purifying selection. It is feasible that weak or sporadic positive selection may occur in this environment with strong purifying selection when lifestyle changes result in increased energy demands or reduced oxygen availability. As a result, phylogenetic reconstruction could take advantage of all PCGs. Besides, the model of evolution among 13 PCGs was mostly in line with the previous literature (Yan et al. 2021b).

Nucleotide diversity among the 13 PCGs is shown in Figure 5. The four PCGs with marked variation were cox1 (Pi = 0.39), nad3 (Pi = 0.354), nad1 (Pi = 0.333), and nad5 (Pi = 0.306), while cox2 (Pi = 0.098), nad4 (Pi = 0.099), cox3 (Pi = 0.106), and cytb (Pi = 0.109) exhibited relatively low Pi values, demonstrating that they are the most conserved genes among the 13 PCGs. Genetic distance analyses also show an analogous tendency (Fig. 5). The mean value of genetic distances within 29 mitogenomes shows that cox1 (mean value = 0.842), nad3 (0.704) and nad1 (0.693) have experienced a comparatively rapid evolution. Inversely, cox2 (0.106), nad4 (0.106) and cox3 (0.115) with lower measured distances are evolving relative slowly.

Overall, we find that the cox1 gene evolves at a considerably higher rate and under comparatively relaxed purifying selection among anthomyiids, manifesting that nad1 gene could be a suitable candidate marker for clarifying the phylogenetic relationships among taxa with morphological traits that are difficult to interpret.

Phylogenetic analyses were carried out on mitogenomes from 29 anthomyiids and eight outgroups, including our 18 newly sequenced Anthomyiidae mitogenomes. Both ML and BI methods were conducted on 13 PCGs and produced completely resolved trees with identical topologies and with most branches receiving strong support. The muscoids were confirmed as a non-monophyletic group or grade, with Scathophagidae plus Anthomyiidae placed as sister to the clade Oestroidea ((Sarcophagidae + Calliphoridae) + Tachinidae), congruent with previous studies (Kutty et al. 2010; Yan et al. 2021a). While the number of clusters within Calyptratae appear to be reliably established, there are still many partial subordinate taxa with weak support or unstable nodes. It is widely known that limited taxon sampling can lead to phylogenetically biased results, and rapid radiation can make resolution of relationships more difficult to resolve (Wiegmann et al. 2011), and this may explain the challenges posed in resolving the muscoid radiation. In our trees, Anthomyiidae was recovered as monophyletic and placed as sister group to the Scathophagidae (BP = 86, PP = 1.00; Fig. 6), in agreement with a recent molecular phylogenetic analyses (Gomes et al. 2021). This finding contrasts with a phylogenomic study using nuclear markers in which Scathophagidae was nested within Anthomyiidae (Buenaventura et al. 2020). In all of these studies, including our own, the Scathophagidae have been represented by limited taxon sampling.

Anthomyiidae and Scathophagidae are clearly very closely related, with various forms of evidence supporting either a sister-group relationship or placement of the latter family within a more broadly defined Anthomyiidae. They are morphologically similar, with both families possessing a long anal vein, usually reaching wing edge at least as a fold (Buck et al. 2009). Reliance on male genitalia to support taxonomy and classification of anthomyiids has made classification of the group technically challenging. Several researchers have provided morphological evidence for the monophyly of Anthomyiidae (Michelsen, 1991; Xue and Chao, 1998), while an increasing number of molecular phylogenetic analyses have found Anthomyiidae to be paraphyletic, containing Scathophagidae (Kutty et al. 2010, 2019). In the Kutty et al. (2010) analysis, the subfamily Anthomyiinae was found to be paraphyletic, comprising 20 species belonging to 10 genera. In addition, Myopininae and Pegomyinae were recovered as monophyletic sister taxa and placed as the sister group of Anthomyiinae, excluding Hydrophoria. Hydrophoria formed a distinct early diverging branch, sister to all other anthomyiids, but with low support (BP = 47, PP = 0.95). Relationships among anthomyiid genera, for example (((Eustalomyia + Leucophora) + Hyporites) + Delia), were strongly supported and as the sister group to the clade ((Botanophila+Hylemyza)+Fucellia) (Kutty et al. 2010). More recently, Gomes et al. (2021) also recovered Botanophila as a sister group of the genus Hylemyza. The specimen labelled Botanophila sp. was difficult to determine morphologically, and here it grouped as sister group to B. fugax. Multiple Pegomya species were polyphyletic, and Emmesomyia oriens emerged as sister to P. flaviprecoxa, with robust support (BP = 100, PP = 1.00). Griffiths (1982) proposed that the condition of a bilobate pregonite is a synapomorphy uniting Pegomya, Emmesomyia and Eutrichota, and this was supported by molecular data (Kutty et al. 2010). It is still noteworthy that several anthomyiid taxa, including Anthomyia, Botanophila, and Delia, are extremely diverse and presumably paraphyletic (Kutty et al. 2008, 2010). Nevertheless, only a single species in each of these genera was used in these studies, thus increased sampling will be necessary to adequately resolve relationships among diverse Anthomyiidae genera.