Adult specimens of Elephantomyia (Elephantomyia) inulta Alexander, 1938 were collected from Motuo, Linzhi, Tibet, China, and adult specimens of Helius (Helius) pluto Alexander, 1932 were collected from Mount Daming, Nanning, Guangxi, China. Specimens were identified based on Zhang et al. (2015b) and Alexander (1932). All specimens were preserved in 96% ethanol at –20 °C for long-term storage at the China Agricultural University, Beijing, China. For each species, genomic DNA was extracted from thoracic muscle tissue using the DNeasy Blood and Tissue kit (Qiagen) according to the manufacturers protocol.

An Illumina TruSeq library was prepared with 450 bp average insert size and sequenced on the Illumina Hiseq 2500 platform with 250 bp paired-end reads. The genomes of the two species were sequenced on one lane. About 4 Gb of clean data was obtained from the library after trimming using Trimmomatic (Bolger et al. 2014). De novo assemblies of high-quality reads were conducted using IDBA-UD 1.1.1 (Peng et al. 2012), with minimum and maximum k values of 80 bp and 240 bp, respectively. Fragments of COI near the 5’-terminus (~610 bp) were amplified for each species by polymerase chain reaction (PCR) with primers LCO1490 (5’-GGTCAACA­AATCATAAAGATATTGG-3’ forward) and HCO2198 (5’-TAAACTTCAGGGTGACCAAAAAATCA-3’ reverse) (Folmer et al. 1994). The PCR was performed in a 25 μL volume containing 12.5 μL Taq PCR Master Mix, 1.0 μL of DNA extract, 1.0 μL primer LCO1490, 1.0 μL Primer HCO2198 and 9.5 μL ddH2O. The cycling profile was 94 °C for 4 min, 30 cycles of 94 °C for 30 sec, 45 °C for 30 sec, 72 °C for 1 min, and a final extension period of 72 °C for 10 min. Successful PCR products were purified and sequenced by Sangon Biotech (Shanghai, China). The COI fragments served as bait references to identify best-fit mt contigs under BLAST searches (Altschul et al. 1990) with minimum similarity 98%. For checking assembly accuracy, clean reads were mapped onto the obtained mt contigs using Geneious 9.0.2 (Kearse et al. 2012).

The protein-coding, rRNA and tRNA genes were identified using the MITOS2 Webserver (http://mitos2.bioinf.uni-leipzig.de/index.py). Protein-coding genes that could not be predicted by that program were annotated by alignment with the homologous genes reported in other crane flies. Nucleotide composition of mt genomes and PCG codon usage was analyzed in MEGA 7.0 (Kumar et al. 2016). AT-skew [(A–T)/(A + T)] and GC-skew [(G–C)/(G + C)] were used to measure nucleotide compositional differences between genes (Perna and Kocher 1995). DnaSP 5.0 (Librado and Rozas 2009) was used to calculate synonymous (Ks) and non-synonymous (Ka) substitution rates, and Ka/Ks ratio (Hurst 2002) was calculated manually.

A total of 31 mt genomes were used for molecular analyses, including the two newly sequenced mt genomes of Elephantomyiini, 26 complete or nearly complete mt genomes of Tipuloidea available in GenBank, and three mt genomes of Trichoceridae used as outgroup (Table 1). The mt genome of Tipula (Pterelachisus) varipennis Meigen, 1818 in Tipuloidea was not included because it lacks both rRNA genes.

The protein-coding and RNA genes were aligned individually with the MAFFT 7.0 online server with the algorithm G-INS-i strategy (Katoh and Standley 2013). Individual gene alignments were checked manually in MEGA 7.0 after removing poorly aligned sites using GBlocks 0.91b (Talavera and Castresana 2007). Some studies suggest that RNA genes evolve too quickly or are too difficult to align to provide phylogenetic signal at deep nodes, but others suggest that the addition of these genes in mt genome phylogenies can improve the resolution and support (Cameron et al. 2007). Alternatively, the third codon positions of PCG have been shown to contain high AT content and accelerated rates of evolutionary change, and so, third codon positions are removed in most phylogenetic studies (Cameron et al. 2009; Mao et al. 2012). However, third position removal also drastically lowers the phylogenetic information content provided by these sites, and so comparison among trees inferred both with and without these 3rd sites included is required to fully evaluate their impact on relationships (Cameron et al. 2007; Fenn et al. 2008; Cameron 2014). Therefore, alignments of individual genes were concatenated using SequenceMatrix 1.8 (Vaidya et al. 2011) to generate the following four datasets: (1) PCGRNA matrix, including all three codon positions of 13 PCGs, large rRNA (lrRNA) and 18 tRNA genes (only these RNA genes are available for all species) (11,353 bp); (2) PCG matrix, including all three codon positions of 13 PCGs (9,444 bp); (3) PCG12RNA matrix, including the first and second codon positions of 13 PCGs, lrRNA and 18 tRNA genes (8,205 bp); and (4) PCG12 matrix, including the first and second codon positions of 13 PCGs (6,296 bp).

Heterogeneous sequence divergence can lead to strong biases in tree reconstructions, such as long branch effects or the misplacement of rogue taxa (Kück et al. 2014), in particular when the taxon sampling is poor, or the outgroup is distant (Felsenstein 1978; Philippe and Laurent 1998). AliGROOVE 1.07, a tool that can detect strongly derived sequences, was used to offer the possibility to exclude taxa or gene partitions (Kück et al. 2014). Phylogenetic trees were inferred for each dataset using Bayesian inference (BI) and maximum likelihood (ML) methods. Partitioning schemes and the best-fit substitution models were determined in PartitionFinder 2.1.1 with the BIC criterion and greedy-search algorithm (Lanfear et al. 2017). The BI analysis was performed in MrBayes 3.2.7 (Huelsenbeck and Ronquist 2001) with the default settings and four independent runs of 5–10 million generations with sampling every 1000 generations; after the average standard deviation of split frequencies fell below 0.01, the initial 25% of samples were discarded as burn-in. The ML analysis was conducted with RAxML 8.2.4 (Stamatakis 2014) using the GTRGAMMAI model; the best ML tree was calculated with branch support estimated from 1000 bootstrap replicates. The program TreePuzzle 5.3 (Strimmer and von Haeseler 1997; Schmidt et al. 2002) was used for four-cluster likelihood mapping (FcLM) analysis to evaluate single topological splits.

The mt genomes of two crane fly species in the tribe Elephantomyiini, E. (E.) inulta and H. (H.) pluto, are sequenced and analyzed for the first time. The nearly complete mt genomes of E. (E.) inulta (GenBank accession no. OP556661) and H. (H.) pluto (GenBank accession no. OP556662) are 14,551 bp and 14,358 bp in length, respectively. The control regions and short stretches on either side of the control regions are not obtained for either species. In the mt genome of E. (E.) inulta, 35 genes are detected (tRNA^Ile, tRNA^Gln and partial small rRNA (­srRNA) are not detected), while in the mt genome of H. (H.) pluto, 34 genes are detected (tRNA^Ile, tRNA^Gln, ­tRNA^Met, partial srRNA and partial ND2 are not detected) (Fig. 3; Table 2). The composition and arrangement in both two mt genomes are identical to the presumed ancestral insect mt genome (Boore 1999). The organization of both mt genomes are generally compact. Intergenic spacers in E. (E.) inulta are 14 in number and generally less than 18 bp, with the largest between tRNA^Ser(UCN) and ND1, while the intergenic spacers in H. (H.) pluto are 19 in number and generally less than 22 bp, with the largest between tRNA^Cys and tRNA^Tyr. Both mt genomes also have overlapping genes but no genes overlapped by more than 8 bp (Table 2).

Figure 3. 

Gene maps of the mitochondrial genomes of two Elephantomyiini species sequenced in this study. The circular maps were drawn with OGDRAW (https://chlorobox.mpimp-golm.mpg.de/OGDraw.html). The transcriptional direction is indicated by arrows.

The mt genomes of both Elephantomyiini species are biased to high A+T% across their four major genome partitions (i.e. PCGs, tRNA genes, lrRNA gene and srRNA gene). The AT contents of whole mt genome, PCGs, tRNA genes and lrRNA gene in E. (E.) inulta (76.4%, 75.2%, 78.8% and 81.5%) are lower than those in H. (H.) pluto (76.8%, 75.8%, 79.4% and 82.1%), but the AT content of the srRNA gene in E. (E.) inulta (78.9%) is higher than that in H. (H.) pluto (76.3%). Both species show slightly positive AT-skew (0.01, 0.02) and negative GC-skew (–0.18, –0.21) for the whole mt genome, but show negative AT-skew (–0.16, –0.16; –0.04, –0.05) and positive GC-skew (0.04, 0.03; 0.33, 0.33) for PCGs and the lrRNA gene. For tRNAs, both species show insignificant or no AT-skew (0.01, 0.00) and positive GC-skew (0.11, 0.14). For the srRNA gene, E. (E.) inulta shows positive AT (0.02) and GC-skews (0.32), while H. (H.) pluto shows negative AT-skew (–0.03) and positive GC-skew (0.27) (Table 3). For both Elephantomyiini species, the whole mt genome and the four major partitions all have the same trends in AT content, AT and GC-skews consistent with the common nucleotide composition of mt genomes of Tipuloidea (Zhang et al. 2016).

Each of the two newly sequenced mt genomes has 13 PCGs, of which COI, COII, COIII, CytB, ATP6, ATP8, ND2, ND3 and ND6 are coded on the majority strand, and ND4, ND4L, ND5 and ND1 are coded on the minority strand (Fig. 3; Table 2). A majority of PCGs in both mt genomes show the typical ATN start codons (ATT/ATG); TCG is the start codon for COI in both species and TTG is the start codon for ND1 in H. (H.) pluto. The majority of PCGs also show the typical TAR (TAA/TAG) stop codons, while the partial stop codon T for COII is found in both species and for ND5 in E. (E.) inulta (Table 2).

The total number of codons of mt genomes are 3,733 in E. (E.) inulta, and 3,700 in H. (H.) pluto but with incomplete ND2 (Tables S1, S2). Codon usage values are described by relative synonymous codon usage (RSCU), which reflects how often each codon is used relative to the expected number in the absence of usage bias. All RSCU values for each amino acid are similar between the two mt genomes, with Leu (UUR) and Ser (UCN) being the two most frequently used amino acids, and Leu (CUN), Met and Trp being the least. The most frequently used codon in each amino acids solely comprises A or T, reflecting the high AT content of PCGs (Fig. S1). These phenomena were also recorded from other mt genomes of lower Diptera (Zhang et al. 2022).

To further investigate evolutionary patterns across the PCGs, the ratio of Ka (rates of nonsynonymous mutations)/Ks (rates of synonymous mutations) are calculated for each (Fig. S2). The Ka/Ks values for all 13 PCGs are lower than 1 (<0.70), implying purifying selection on all these genes. The Ka/Ks ratio of ND2 is obviously higher than other PCGs, which indicates that ND2 has a relatively higher evolutionary rate. In contrast, COI has the lowest Ka/Ks ratio, indicating that this gene has been subjected to the highest purifying selection.

Twenty tRNA genes are detected in E. (E.) inulta and 19 tRNAs in H. (H.) pluto. The tRNA genes lengths range from 64 bp to 72 bp (Table 2). Most tRNAs can be folded into the typical cloverleaf structure, except for tRNA^Ser(AGN) whose dihydrouridine (DHU) arm is replaced by a simple loop (Fig. S3). It is common for tRNA^Ser(AGN) to lack the DHU arm in insect mt genomes (Zhang et al. 2016; Zhang et al. 2018; Zhu et al. 2018; Ren et al. 2019; Su and Liang 2019; Wang and Huang 2019; Li et al. 2021; Mo et al. 2022; Zhang et al. 2022). Nucleotide substitutions of tRNAs between E. (E.) inulta and H. (H.) pluto range from three to 23: tRNA^Asp has the least variation with three substitutions, while tRNA^Ala has the most variation with 23 substitutions and indels.

As in the ancestral insect (Zhang et al. 2016, 2022; Chen et al. 2017; Liu et al. 2017; Feng et al. 2018; Zhang et al. 2018; Zhu et al. 2018; Ren et al. 2019; Su and Liang 2019; Wang and Huang 2019; Zhao et al. 2019; Tang et al. 2020; Li et al. 2021; Wang et al. 2021; Shi et al. 2021; Sun et al. 2021; Zheng et al. 2021; Mo et al. 2022), the lrRNA gene in each Elephantomyiini species is located between tRNA^Leu(CUN) gene and tRNA^Val gene, while the srRNA gene is located between tRNA^Val gene and the control region (not obtained in this study) (Fig. 3). The lengths of lrRNAs are 1,326 bp in E. (E.) inulta and 1,321 bp in H. (H.) pluto. The assembled srRNA gene in both species are incomplete, and the obtained sequence lengths are 615 bp and 558 bp respectively (Table 2).

AliGROOVE analysis indicates that Chionea crassipes gracilistyla Alexander, 1936 has the strongest heterogeneity relative to other Tipuloidea species in all four datasets (Fig. S4), which may cause bias tree reconstructions and node support in phylogenetic analysis (Kück et al. 2014). An effective solution is to avoid using this species (Soltis et al. 2004). Therefore, C. crassipes gracilistyla was excluded when constructing phylogenetic trees. Twenty-eight representatives from all four families of Tipuloidea and three representatives from Trichoceridae were included in the phylogenetic analysis. Eight phylogenetic trees inferred from the four datasets under BI and ML methods were finally constructed (Figs 4, S5–S10), resulting in four hypotheses of relationships among major groups of Tipuloidea (Fig. 5).

Figure 4. 

Phylogenetic trees of Tipuloidea inferred from the datasets A PCG12RNA and B PCG under BI method. Numbers at the nodes are posterior probabilities. The family Trichoceridae was set as the outgroup.

Figure 5. 

Four hypotheses for the relationships among major groups of Tipuloidea in this study. Multiple sampling of different species from a single family/subfamily/tribe are collapsed into triangles. Numbers at the nodes are posterior probabilities/bootstrap values, NS = not support.

In all BI and ML trees, Pediciidae is sister to all other Tipuloidea, and a sister relationship between Cylindrotomidae and Tipulidae is strongly supported. These arrangements are consistent with the phylogeny by Ribeiro (2008) based on 88 morphological characters, Petersen et al. (2010) based on combined morphological characters and two nuclear genes, Zhang et al. (2016) based on mt genomes and Kang et al. (2017) based on transcriptomes.

Limoniidae is not supported as monophyletic clade in any phylogenetic trees. Symplecta (Symplecta) hybrida (Meigen, 1804) (Chioneinae) is sister to all non-pediciid crane flies in trees inferred from the PCG12RNA and PCG12 datasets under BI and ML methods (96% PP, 45% BV; 64% PP, 42% BV) (Figs 4A, S5, S7, S8), while in the BI tree inferred from the PCG dataset (Fig. 4B), Symplecta is sister to a clade of Cylindrotomidae + Tipulidae (93% PP). In the both BI and ML trees inferred from the PCGRNA dataset (Figs S9, S10), Symplecta + the Limnophilinae species Pseudolimnophila (Pseudolimnophila) brunneinota Alexander, 1933 is weakly supported (55% PP, 60% BV), as sister to all non-pediciid crane flies (54% PP, 33% BV). In the ML tree inferred from the PCG dataset (Fig. S6), Symplecta is sister to a clade of Limnophilinae + (Cylindrotomidae + Tipulidae) with a very low bootstrap value (19% BV). However, our current taxon sampling for Chioneinae is far from extensive, and further detailed studies with more taxa are needed before the monophyly of Chioneinae can be confidently defined.

Limoniinae (including Elephantomyiini) + Epiphragma (Epiphragma) mediale Mao and Yang, 2009 (Limnophilinae) forms a clade in all phylogenetic trees (100% or 99% PPs for all BI trees; 65%, 56%, 43% and 37% BVs for ML trees). Elephantomyiini (Elephantomyia + Helius) (100% PP for all BI trees; 94%, 92%, 92% and 89% BV for ML trees) and Limoniinae (100% PP/BV for all trees) are two well-supported clades, which to some extent supports the suggestion of Petersen et al. (2010) to treat the monophyletic Elephantomyiini as a subfamily on the basis of morphological and molecular data. Although Epiphragma consistly shows a distant relationship to other Limnophilinae in our study and was also treated as a subfamily in Petersen et al. 2010, the taxonomic status of the group it represents needs further study.

Limnophilinae is a controversial group with respect to both its monophyly and relationships with other Limoniidae. The main clade of Limnophilinae (including four species) is sister to the clade containing Elephantomyiini, Epiphragma and Limoniinae in the trees inferred from the PCG12RNA dataset under BI and ML methods (95% PP, 27% BV) (Figs 4A, S5), while in the remaining trees, a clade of Limnophilinae (including three or four species) shows a closer relationship to Cylindrotomidae + Tipulidae.

Topologies I and II (Fig. 5) based on the PCG12RNA and PCG datasets under BI method have high PPs (Fig. 4), but are contradictory. FcLM analysis is often used to evaluate single topological splits (Trautwein et al. 2010; Misof et al. 2014; Peters et al 2014; Winterton and Ware 2015; Kang et al. 2017; Narayanan et al. 2018, 2019; Vasilikopoulos et al. 2019; Zhang et al. 2019; Karmeinski et al. 2021; Wang et al. 2022). Two questions that reflect the main differences between topologies I and II were used for FcLM testing to further evaluate these two topologies: 1) Is Symplecta sister to all non-pediciid crane flies, or to Cylindrotomidae + Tipulidae? 2) Is Limnophilinae part of Limoniinae, or does Limnophilinae have a closer relationship with Cylindrotomidae + Tipulidae? (Table 4). For these tests, species in four datasets were grouped into four clusters representing alternative resolutions of these topologies.

Our FcLM analysis shows a support for the sister-group relationship between Symplecta and all non-pediciid crane flies (51.0%/43.9%/81.4%/81.3%) (Fig. S11). FcLM results for the placement of Symplecta are concordant with topology I (Fig. 5). However, the status of Limoniinae (except Chioneinae) is not well resolved by FcLM analysis (Fig. S12): Limnophilinae is supported as part of Limoniinae (as shown in topology I) with the support rate of 40.8%/31.4%/36%/18.1% whereas Limnophilinae is supported as sister to Cylindrotomidae + Tipulidae (as shown in topology II) with the support rate of 30.5%/31.8%/40.4%/49.7%.