All the specimens of the three species sequenced in this study were collected and identified by authors and the collecting information is summarized in Suppl. material 1. After collection, each specimen was immediately preserved in 95% ethanol, and later stored at -20 °C in the laboratory. The voucher specimen of T. (P.) cinereocincta mesacantha was deposited in the Entomological Museum of Shenyang Agricultural University (SYAU), while those of T. (P.) legalis and T. (P.) yasumatsuana were deposited in the Entomological Museum of China Agricultural University (CAU). Genomic DNA was extracted from the thoracic tissue of each specimen using the QIAamp DNA Blood Mini Kit (Qiagen, Germany). The DNA concentration was quantified using an Agilent 5400 instrument.

For DNA library preparation, the NEB Next® Ultra™ DNA Library Prep Kit was utilized, and paired-end sequencing was conducted on an Illumina NovaSeq 6000 platform, generating raw data with an insert size of 350 bp and a read length of 150 bp. Approximately 4 Gb of raw sequenced data was obtained. Novogene Biotechnology Company (Beijing, China) conducted the aforementioned processes. The paired raw reads for the whole mitogenome of T. (P.) varipennis was downloaded from NCBI under the accession number SRR11469981.

The mitochondrial genomes of all species were assembled using IDBA-UD 1.1.3 (Peng et al. 2012), and circularization of resulting linear contigs was verified using the python script in MitoZ 2.3 software (Meng et al. 2019). Circular mitochondrial genomes were then submitted to the MITOS2 web service (Bernt et al. 2013) for annotation. The secondary structure of tRNA was determined using both the MITOS2 web service (Bernt et al. 2013) and the tRNAscan-SE 2.0 (Lowe and Chan 2016) web service. Annotated mitochondrial genomes underwent a comparative analysis with closely related species in Geneious 9.0.2, and manual corrections were applied. Subsequent analyses were conducted after error exclusion.

Gene maps of the mitochondrial genomes sequenced of four T. (Pterelachisus) species were generated using the Proksee web service (Grant et al. 2023). Basal composition and amino acid usage were calculated using PhyloSuite 1.2.3 (Zhang et al. 2020), with AT skew defined as [A – T] / [A + T] and GC skew defined as [G – C] / [G + C] (Perna and Kocher 1995). Ka and Ks values, along with nucleotide diversity (Pi), were obtained using DnaSP 6.12.03 (Rozas et al. 2017). Relative synonymous codon usage (RSCU) data were also acquired through PhyloSuite (Zhang et al. 2020), and Rscript was employed for graphical representation. Repeat segments in the control region (CR) were identified using the Tandem Repeats Finder 4.09 (Benson 1999).

Accurately annotated mitochondrial genomes, along with sequencing data, were deposited in the NCBI database under the BioProject PRJNA1067446.

A total of 31 complete mitochondrial genomes were used for phylogenetic analysis in this study (Table 1). Trichocera bimacula Walker, 1848 and Paracladura trichoptera (Osten Sacken, 1877), members of the family Trichoceridae, were designated as outgroups, serving as the root of the phylogenetic tree. Twenty-nine species within four families of Tipuloidea were contained in the ingroups, which respectively include one species of Pediciidae, ten species within ten genera of Limoniidae, one species of Cylindrotomidae, and 17 species within 12 (sub)genera of Tipulidae. All data preprocessing was carried out using PhyloSuite 1.2.3. Mitochondrial genomes served as the basis for constructing four concatenated datasets: 1) 13PCG, including all three codon positions of 13 PCGs. 2) 13PCG + rRNA, including all three codon positions of 13 PCGs and two Ribosomal RNA genes. 3) 13PCG12, including the first and second codon positions of 13 PCGs. And 4) AA, including all amino acid of 13 PCGs. Prior to concatenation, all data underwent alignment using MAFFT (Katoh and Standley 2013), followed by manual correction in MEGA7 (Kumar et al. 2016) to eliminate gap regions. Model selection for optimal models was performed using PartitionFinder 2.1.1 (Lanfear et al. 2016).

AliGROOVE 1.08 (Kück et al. 2014) was used to offer the possibility to exclude taxa or gene partitions. Phylogenetic analysis for Maximum Likelihood (ML) trees utilized RAxML 8.2.12 (Stamatakis 2014) with specific parameters set as -m GTRGAMMA -x 1234 -p 12345 -# 1000. Bayesian analysis was conducted using MrBayes 2.3 (Ronquist et al. 2012) for 2,000,000 generations with the default settings. The resulting phylogenetic tree was visualized and enhanced for presentation using Figtree 1.4.4 and Adobe Photoshop 2022.

The complete mitochondrial genomes of all four T. (Pterelachisus) species comprise 13 protein-coding genes, 22 transfer RNA genes, two ribosomal RNA genes, and one non-coding region (A + T-rich control region) (Table 2; Fig. 1). These genes exhibit a ring structure and the tandem arrangement is consistent with previously published mitochondrial whole genome gene arrangements in species of Tipulidae. The total length of the four mitochondrial genomes ranges from 15,000 to 16,000 base pairs. Specifically, T. (P.) cinereocincta mesacantha, T. (P.) legalis, T. (P.) varipennis, and T. (P.) yasumatsuana have lengths of 15,907 bp, 15,625 bp, 15,772 bp, and 15,735 bp, respectively (Table 3). All mitochondrial genomes are notably AT-rich, with A + T base contents for T. (P.) cinereocincta mesacantha, T. (P.) legalis, T. (P.) varipennis, and T. (P.) yasumatsuana at 76.7%, 75.0%, 77.8%, and 75.4%, respectively (Table 3).

Figure 1. 

Gene maps of the mitochondrial genomes of the four T. (Pterelachisus) species involved in this study. The transcriptional direction is indicated by arrows.

The mitochondrial genomes of the four species share similar, but not identical, intergenic regions and overlaps. The longest intergenic regions, found between trnE and trnF genes, measure 32 bp, 27 bp, 32 bp, and 27 bp for T. (P.) cinereocincta mesacantha, T. (P.) legalis, T. (P.) varipennis, and T. (P.) yasumatsuana, respectively. The longest overlaps, located between trnW and trnC genes, are consistent across all species at a length of 8 bp.

All four mitochondrial genomes harbor 13 protein-coding genes, including COX1, COX2, COX3, CYTB, ATP6, ATP8, ND2, ND3, and ND6 on the majority strand, and ND4, ND4L, ND5, and ND1 on the minority strand (Fig. 1; Table 2). All species exhibit a pronounced AT richness, with A + T base content for T. (P.) cinereocincta mesacantha, T. (P.) legalis, T. (P.) varipennis, and T. (P.) yasumatsuana at 74.6%, 72.6%, 76.0%, and 73.0%, respectively. The AT richness is especially evident in third codon positions, all exceeding 80.0%, with T. (P.) varipennis having the highest value at 92.4%. The first and second codon positions have lower AT skewness values, all below 70.0%. The most frequently encoded amino acids in these four T. (Pterelachisus) mitogenomes are Ser2, Leu2, Val, Gly, Pro, Thr, Arg, and Ala, with the highest Relative Synonymous Codon Usage (RSCU) values (Fig. 2). The most common codons are UUA, AUU, UUU, and AUA, and the majority of codons are composed solely of A or T, reflecting the high AT content of protein-coding genes (PCGs).

Figure 2. 

Relative synonymous codon usage (RSCU) in the mitogenomes of the four T. (Pterelachisus) species. Codes as follows: A: Ala; C: Cys; D: Asp; E: Glu; F: Phe; G: Gly; H: His; I: Ile; K: Lys; L: Leu; M: Met; N: Asn; P: Pro; Q: Gln; R: Arg; S: Ser; T: Thr; V: Val; W: Try; Y: Tyr.

For most PCGs, typical ATN start codons (ATT / ATG) are observed in both mitochondrial genomes, except for TCG in COX1 genes. Stop codons for most PCGs are T + tRNA, while CYTB has a stop codon TAG (Table 2). The sliding window analysis reveals variable nucleotide diversity (Pi) among the 13 PCGs in the four mitochondrial genomes, with ND2 exhibiting the highest Pi (0.308), followed by ATP8 (0.276) and ND6 (0.257). ND5 shows the lowest Pi (0.151) (Fig. 3A). Further examination of the Ka / Ks ratio for each PCG indicates values less than 1, suggesting purification selection. ND2 has a notably higher Ka / Ks ratio, indicating a higher evolutionary rate, while COX1 undergoes the highest purification selection. The Ka / Ks ratio of ND6 varies significantly among the four species, with T. (P.) varipennis having significantly higher values than the other species (Fig. 3B, Suppl. material 2).

Figure 3. 

A The nucleotide diversity (Pi) of 13 protein-coding genes (PCGs) in four T. (Pterelachisus) species mitogenomes determined via sliding window analysis (sliding window: 100 bp; step size: 25 bp); the Pi value of each gene is shown under the gene name B evolutionary rates (ratios of Ka/Ks) of mitochondrial protein-coding genes of the four T. (Pterelachisus) species.

All mitochondrial genomes encompass 22 tRNA genes, each capable of forming cloverleaf structures, with the exception of trnS1 (AGC), which has a dihydrouridine (DHU) arm forming a loop (Fig. 4). The length of the 22 tRNA genes ranges from 62 to 72 bp across the four mitochondrial genomes. The shortest trnC (GCA) genes are found in all species, with a length of 62 bp, except for trnC in T. (P.) varipennis, which is 63 bp. The longest tRNA genes are trnV (CAC) genes.

Figure 4. 

Secondary structures of tRNAs of T. (P.) varipennis. All tRNAs are labeled with the abbreviations of their corresponding amino acids. The variable sites are indicated with the green coloration for T. (P.) cinereocincta mesacantha, blue coloration for T. (P.) legalis and pink coloration for T. (P.) yasumatsuana, respectively. A blank represents a missing base site.

The tRNA genes in all four species exhibit significant AT richness, with A + T base content for T. (P.) cinereocincta mesacantha, T. (P.) legalis, T. (P.) varipennis, and T. (P.) yasumatsuana at 76.5%, 77.2%, 77.4%, and 76.8%, respectively (Table 3).

All four mitochondrial genomes feature two ribosomal RNA genes, rrnL and rrnS, separated by trnV. The rrnL of all four species (1,319 bp –1,323 bp) is notably longer than the rrnS (782 bp –783 bp). The rRNA is significantly AT-rich in all species, with A + T base content for T. (P.) cinereocincta mesacantha, T. (P.) legalis, T. (P.) varipennis, and T. (P.) yasumatsuana at 80.2%, 79.1%, 80.4%, and 79.4%, respectively (Table 3).

The control region for all four species is situated between rrnS and trnI genes, with lengths ranging from 800 bp to 1,100 bp. T. (P.) cinereocincta mesacantha has the longest control region at 1,806 bp, while T. (P.) legalis has the shortest at 805 bp. The control regions of all four species exhibit significant AT richness, with T. (P.) cinereocincta mesacantha, T. (P.) legalis, T. (P.) varipennis, and T. (P.) yasumatsuana having A + T base content of 90.7%, 91.5%, 91.4%, and 91.4%, respectively (Table 3).

The control regions for all four species were analyzed using the Tandem Repeat Finder, revealing two or three tandem repeats of varying lengths (Fig. 5). T. (P.) cinereocincta mesacantha and T. (P.) varipennis have two sets of tandem repeats, while T. (P.) legalis and T. (P.) yasumatsuana have three sets. Tandem repeats between different species exhibit no obvious common features, displaying unique structural and evolutionary characteristics.

Figure 5. 

Control region structure of the four T. (Pterelachisus) species. Orange, pink, and red coloration represent tandem repeats.

Both Bayesian inference (BI) and Maximum Likelihood (ML) trees were reconstructed using four concatenated datasets (13PCG, 13PCG12, 13PCG + rRNA, and AA) of 31 mitochondrial genomes (Fig. 6, Suppl. material 3: figs S1–S6). The topologies of these trees show notable differences. Heterogeneity in pairwise sequence differences was examined and the results have shown the AA dataset with significantly lower heterogeneity compared to the other datasets (Fig. 7), which may be a major factor influencing the phylogenetic results.

Figure 6. 

Phylogenetic trees of the selected species of Tipuloidea inferred from the datasets PCG under A ML and B BI methods. Numbers at the nodes are bootstrap values (ML tree) or posterior probabilities (BI tree). The two species of family Trichoceridae were set as the outgroups.

Figure 7. 

AliGROOVE analysis for four datasets. The mean similarity score between sequences is represented by a colored square, based on AliGROOVE scores ranging from -1, indicating a large difference in sequence composition from the remainder of the dataset (red coloration), to +1, indicating similarity to all other comparisons (blue coloration).

The four T. (Pterelachisus) species involved in this study are divided into two stable lineages in each phylogenetic tree above: T. (P.) cinereocincta mesacantha and T. (P.) varipennis form a sister group, while T. (P.) legalis and T. (P.) yasumatsuana form another one. Furthermore, almost all the trees, except those based on the AA dataset (Fig. 7, Suppl. material 3: figs S1, S2, S4, S5), have shown that the four T. (Pterelachisus) species compose a monophyletic lineage, but with variable support values among different trees, whereas both the ML and BI trees based on the AA dataset suggest T. (Pterelachisus) is a paraphyletic group (Suppl. material 3: figs S3, S6). Since the samples for the large subgenus T. (Pterelachisus) used in this study are far from sufficient, the monophyly of T. (Pterelachisus) needs further study with more data.

The above “Vestiplex-Lunatipula” group of the Tipula subgenera are tentatively supported to be a monophyletic lineage by the phylogenetic results based on the 13PCG dataset (Fig. 6), but then the subgenus T. (Formotipula), unexpectedly, should be included in this group, in addition to T. (Vestiplex), T. (Lunatipula) and T. (Pterelachisus). These arguments largely agree with the phylogenetic results based on the 13PCG + rRNA dataset (Suppl. material 3: fig. S5) in the present study, as well as the previous research results of Gao et al. (2023), only with different topologies. On the contrary, the phylogenetic trees based on the AA (Suppl. material 3: figs S3, S6) and 13PCG12 (Suppl. material 3: figs S1, S4) datasets show that the “Vestiplex-Lunatipula” group is paraphyletic, and T. (Vestiplex) is a sister-group to the remaining Tipulidae. Morphologically, T. (Pterelachisus) and other related subgenera such as T. (Vestiplex), T. (Lunatipula), and other subgenera of Tipula share many similarities, making it difficult to distinguish between them (Savchenko 1964; Gelhaus 1986; 2005). This supports their potential monophyly to some extent. However, it is challenging to explain the close phylogenetic relationship of T. (Formotipula) with these subgenera. Obviously, in-depth research would be required to resolve the questions on its monophyly and relationships within the group. In addition, the genus Tipula is not supported to be a monophyletic lineage by any of the phylogenetic trees.

The monophyly of Tipulidae and the sister relationship between Tipulidae and Cylindrotomidae are strongly supported in all BI and ML trees constructed in this study, which are consistent with the previous phylogenetic studies of Ribeiro (2008), Petersen et al. (2010), Zhang et al. (2016), and Kang et al. (2017; 2023). Tipulidae was divided into three subfamilies, i.e., Ctenophorinae, Dolichopezinae, and Tipulinae (Kertesz 1902; Oosterbroek 2024). In the present phylogenetic study, Dolichopezinae is not included due to a lack of mitogenomic data on the subfamily. Meanwhile, as the only representative of Ctenophorinae, Tanyptera hebeiensis Yang & Yang, 1998 is sister to some members of Tipulinae (i.e., Tipula (Dendrotipula) flavolineata Meigen, 1804, Tipula (Vestiplex) aestiva Savchenko, 1960, or Nephrotoma spp. in different topologies) and then clustered with some other members of Tipulinae, which indicates a para- or polyphyly of the subfamily Tipulinae.

Corroborating previous phylogenetic studies (Ribeiro 2008; Petersen et al. 2010; Zhang et al. 2016; Kang et al. 2023), Limoniidae is confirmed as a non-monophyletic group in this study. Among the four traditional subfamilies of Limoniidae established by Starý (1992), three are involved in the present phylogenetic study, i.e., Chioneinae, Limnophilinae, and Limoniinae. With relatively low supporting values, Symplecta hybrida (Meigen, 1804), the only representative of Chioneinae in this study, is sister to a clade of Cylindrotomidae + Tipulidae in almost all the trees except the BI one inferred from the dataset AA (Suppl. material 3: fig. S6). The traditional subfamily Limnophilinae is not supported as a monophyletic group, because one of its members, Epiphragma mediale Mao & Yang, 2019, has a relatively stable sister relationship with the clade of Limoniidae, instead of with other species of Limnophilinae, which is shown in all the BI trees and most ML trees except the one inferred from the dataset AA (Suppl. material 3: fig. S5). Similar results were also indicated in the previous studies by Kang et al. (2023) and Xu et al. (2023). Furthermore, the family Pediciidae is well supported to be sister to the remaining Tipuloidea, as all the phylogenetic analyses available on the Tipuloidea since Starý (1992).

The authors have declared that no competing interests exist.

No ethical statement was reported.

This research was supported by the National Natural Science Foundation of China (31970444; 31501880) and the Scientific Research Foundation for the Introduced Talent of Shenyang Agricultural University (880415013).

Yan Li and Ding Yang planned and designed the research. Yuetian Gao and Wanxin Cai performed experiments, and Yuetian Gao and Yupeng Li analyzed the data. Yuetian Gao wrote and other authors revised the manuscript. All authors have approved the manuscript for publication and agreed to be accountable for all aspects of the work.

Yuetian Gao https://orcid.org/0000-0003-2966-1745

Wanxin Cai https://orcid.org/0000-0001-6604-4899

Yupeng Li https://orcid.org/0009-0008-5381-2260

Yan Li https://orcid.org/0000-0002-4896-5843

Ding Yang https://orcid.org/0000-0002-7685-3478

All of the data that support the findings of this study are available in the main text or Supplementary Information.