Comparative Analysis of Mitogenomes of Microtendipes (Diptera: Chironomidae) and Its Phylogenetic Implications

1. Introduction

Insect mitochondrial genomes are typically circular molecules ranging in size from 14 to 20 kilobases (kb), encoding a conserved arrangement of 13 protein-coding genes (PCGs), 22 transfer RNAs (tRNAs), and two ribosomal RNAs (rRNAs) in a characteristic order and orientation across most insect orders [1]. These genomes have proven instrumental in diverse fields, including species identification [2,3,4] and population genetics [5,6]. This is due to their advantages, which stem from maternal inheritance, rapid substitution rates, and ease of accessibility [7]. Moreover, investigating mitochondrial genome characteristics, such as nucleotide composition, codon usage, evolutionary rates, and secondary structures of RNA genes, has significantly contributed to a deeper understanding of the evolutionary trajectories of diverse organismal groups [8,9,10,11,12,13]. Within the insect family Chironomidae, mitogenomes have been extensively utilized in studies aimed at elucidating phylogenetic relationships within the family and exploring its evolutionary history [14,15,16,17].

Chironomidae, normally known as non-biting midges, exhibit the greatest diversity of species that have evolved to thrive under extreme abiotic environments. These include species capable of enduring a wide range of temperatures, from extremely low temperatures to extremely high temperatures, as well as those resilient to oxygen scarcity, high levels of salinity, both acidic and alkaline conditions (low and high pH), and even complete desiccation [18,19]. Their remarkable adaptability allows them to inhabit an extensive range of habitats, spanning from the frigid, glacier-covered peaks of the tallest mountains to the depths of freshwater bodies, demonstrating their remarkable versatility and resilience in diverse environments [20,21].

Microtendipes Kieffer, 1915, a genus of the tribe Chironomini in the subfamily Chironominae, exhibits a cosmopolitan distribution across all zoogeographical regions [22]. Its immature stages are predominantly found within littoral and sublittoral sediments of vast aquatic habitats, with a limited number of species also inhabiting flowing water environments. The extensive variability in body colors, encompassing thoracic, leg, and wing pigmentation, poses challenges in defining species solely based on color pattern changes [23]. DNA barcode data are pivotal in supporting the interpretation of pigmentation variations as interspecific differences [24].

Previous studies have employed morphological analysis and molecular markers to investigate species delimitation and phylogeny within Microtendipes [22,23,24], yet certain questions remain unaddressed. Firstly, a complex scenario arises when samples displaying identical color patterns form distinct clades on the evolutionary tree based on COI (cytochrome oxidase subunit I) gene sequences, accompanied by significant genetic distances, rendering species delimitation uncertain. There is debate regarding the appropriateness of classifying deeply divergent specimens, forming paraphyletic groups in COI-based phylogenetic trees, as the same species, as suggested by Song et al. (2024) [24]. This underscores the need for unexplored avenues such as the application of mitogenomes to elucidate phylogenetic relationships. Prior to this investigation, only a single Microtendipes species had been documented with a mitogenome [25], and no comparative analysis of nucleotide composition or evolutionary rates within the genus had been conducted. Secondly, the phylogenetic position of Microtendipes has long been a subject of controversy. Sæther (1977) constructed a phylogenetic tree based on female chironomid characteristics, placing this genus within Chironomini [26]. In the phylogenetic tree, Microtendipes formed a monophyletic group with Paratendipes Kieffer /Nilothauma Kieffer, leading to the suggestion of establishing these genera as a fourth subfamily-level taxonomic unit. However, when Sasa (1989) documented the chironomid fauna of Japan, he placed Microtendipes within the Polypedilum complex (including Nilothuma and Paratendipes), belonging to the Chironomini [27]. When molecular data were introduced, Cranston et al. (2012) supported the grouping of Microtendipes with several related genera (excluding Paratendipes) into a generic group, which formed a sister relationship with Pseudochironomini [28].

In this study, we have presented eight specimens belonging to five species of Microtendipes, specifically M. bimaculatus, M. baishanzuensis, M. wuyiensis, M. robustus, and M. tuberosus. Following this, we conducted a thorough analysis of their mitogenomic characteristics. By incorporating previously documented mitogenomes into our analysis, we reconstructed the phylogenetic evolution of the Chironominae using a comprehensive dataset comprising 25 mitochondrial genomes. To achieve this, we employed both Maximum Likelihood (ML) and Bayesian Inference (BI) methodologies, enabling us to accurately determine the phylogenetic positioning of Microtendipes within the subfamily.

2. Materials and Methods

2.1. Taxon Sampling and Sequencing

All the samples of Microtendipes individuals were collected from Baishanzu National Reserve and Wuyi Mountain National Reserve in the year of 2021 and 2022 except Microtendipes bimaculatus-Clade1 (MZ981734) retrieved from our previous study [25]. Targeted specimens were dissected and mounted into Euparal except the tissues for total genomic DNA extraction (thorax, head and a pair of legs). The extraction procedure followed the Qiagen DNeasy Blood and Tissue kit except for elusion buffer ranging from 100-150 µl according to different body size. The exoskeletons of thorax were cleared and mounted on corresponding voucher. A segment of cox1-5P for each species was amplified with primer pair LCO1490-L and HCO2198-L [29] and sequenced with Sanger sequencing as described by Song et al. (2018) [30], and further confirmed the identification. The genomic DNA was subsequently pooled and sequenced using the Illumina Novaseq 6000 (PE150, Ullumina, San Diego, CA) platform. Multi softwares, Mitoz (v3.4) [31], NOVOPlasty (v4.3.1) [32] and SPAdes (v4.0.0) [33], were employed to assemble the reads, and the results of three assemblers were used to revise each other to obtain accurate mito genomes. The mito genomes were annotated using MITOS2 (v2.1.7) [34] and then visualized using java package CGView (v2.0.3) [35].

2.2. Genome Composition and Codon Usage

The base composition was analyzed utilizing seqkit (v2.3.0) [36]. The codon usage of PCGs were analyzed using EMBOSS (v6.6.0.0) [37] and python module codonw-slim (v1.5.0). The figures about codons of PCGs were generated using R package ggplot2 (v3.5.1) [38].

2.3. Substitution Rate and Phylogenetic Analyses

The mitochondrial genomes of species listed in Table 1, including the eight newly sequenced mitogenomes reported in this study, were selected for phylogenetic reconstruction. Rheocricotopus villiculus (MW373526) was designated as the outgroup. The thirteen protein-coding genes (PCGs) were aligned using MAFFT (version 7.505) [39], with all stop codons retained and indels removed during the process. Subsequent to alignment, all sequences were manually reviewed, corrected, and concatenated for further analyses. Pairwise genetic distances among the 13 PCGs were computed using MEGA11 [40]. Statistical comparisons of genetic distances were conducted using the Friedman test in R, with multiple pairwise comparisons performed using the Wilcoxon rank-sum test, and p-values adjusted using the Benjamini-Hochberg (BH) method to control for false discovery rate. The resulting boxplot was generated using the R package ggplot2.

A total of four datasets were concatenated for phylogenetic analyses, namely: (i) a PCG matrix comprising all thirteen protein-coding genes (PCGs) totaling 11,094 base pairs (bp), (ii) a PCGrRNA matrix incorporating the thirteen PCGs and the two ribosomal RNA (rRNA) genes, amounting to 12,643 bp, (iii) a PCG12 matrix containing the first and second codon positions of the thirteen PCGs, and (iv) a PCG12rRNA matrix (7,396 bp) that includes PCG12 and the two rRNA genes, totaling 8,945 bp. All matrices underwent analysis using maximum likelihood (ML) with RAxML (version 8.2.12) [41] and Bayesian inference (BI) with MrBayes (version 3.2.7a) [42]. PartitionFinder2 (version 2.1.1) [43] was employed to identify the optimal partitioning schemes and substitution models. Both ML and BI analyses were performed on both partitioned and non-partitioned data to construct the phylogenetic trees. For the ML analysis, the "GTRGAMMA" nucleotide substitution model was selected, and a bootstrap analysis with 1000 replicates was conducted. In the BI analysis, four simultaneous Markov chain Monte Carlo (MCMC) runs of 100 million generations each were executed, with trees sampled every 2000 generations. The first 25% of the steps were discarded as burn-in. Stationarity was deemed achieved when the average standard deviation of split frequencies dropped below 0.01. The resulting trees from all analyses were visualized and edited using iTOL (https://itol.embl.de) [44].

3. Results

3.1. Genome Organization

The complete mitochondrial genomes of Microtendipes species have been successfully sequenced. These genomes consist of double-stranded circular molecules, with sizes ranging from 15,756 bp (PP966950) to 16,340 bp (PP966952). Each genome encodes a total of 37 genes, comprising 13 protein-coding genes (PCGs), 2 ribosomal RNA (rRNA) genes, and 22 transfer RNA (tRNA) genes, as detailed in Table 1. The primary source of length variation among these genomes is the control region, which exhibits a size range from 338 bp (PP966951) to 1,266 bp (PP966952). Regarding the coding sequences, minimal variation is observed in the lengths of the PCGs, tRNAs, and rRNAs. Specifically, the PCGs span from 11,106 bp (PP966952) to 11,154 bp (PP966949), the tRNAs range from 1,575 bp (PP966952) to 1,595 bp (PP966949), and the rRNAs vary between 1,488 bp (PP966950) and 1,508 bp (PP966949). Notably, these 37 genes maintain a consistent order across all nine genomes examined. Furthermore, 23 of these genes (including 9 PCGs and 14 tRNAs) are encoded by the majority strand (J-strand), while the remaining genes are encoded by the minority strand (N-strand), as illustrated in Figure 1.

3.2. Nucleotide Composition

The nucleotide compositions of the nine mitochondrial genomes exhibit considerable similarity, as outlined in Table 2. The AT content varies between 79.19% in Microtendipes wuyiensis (PP966951) and 80.77% in Microtendipes robustus (PP966954). Notably, the AT content within the control region is significantly elevated compared to the entire genome, with a range from 94.14% in Microtendipes bimaculatus (PP966953) to 96.42% in Microtendipes robustus (PP966950). Similarly, the AT content in both rRNA and tRNA regions surpasses the overall genome level, with rRNA varying between 81.43% and 82.51%, and tRNA ranging from 82.18% to 82.93%. Conversely, the PCG sequences demonstrate a lower AT content, spanning from 76.07% in Microtendipes wuyiensis to 78.09% in Microtendipes robustus (PP966954). All nine mitogenomes exhibit a positive AT-skew, ranging from 1.08% to 2.66%, and a negative GC-skew, varying between -22.26% and -18.04%. In contrast, the tRNA sequences display a negative AT-skew, ranging from -2.43% to -3.76%, and a positive GC-skew, spanning from 27.47% to 31.67%. The rRNA sequences exhibit positive AT-skew and GC-skew values, with ranges of 3.92% to 5.05% and 12.23% to 16.23%, respectively. For the PCG sequences, a notable negative AT-skew is observed, ranging from -17.31% to -19.16%, whereas the GC-skew is relatively insignificant, varying from -1.57% to 2.38%.

3.3. Protein-Coding Genes and Codon Usage

Among the 13 protein-coding genes (PCGs) in the nine mitochondrial genomes, 12 utilize ATG/ATT as their start codons. An exception is the cox1 gene, which employs ATA as the start codon in Microtendipes wuyiensis (PP966951) and TTG in the other eight mitogenomes. The nad1 gene in Microtendipes wuyiensis terminates with a TAG stop codon, whereas all other genes in the remaining mitochondrial genomes utilize TAA as their stop codon. Tables S1 and S2 provide the amino acid compositions and codon usage tables, respectively. Excluding stop codons, the total number of codons ranges from 3,689 to 3,705. Notably, the four M. bimaculatus mitogenomes contain 3,689, 3,689, 3,694, and 3,694 codons. The most frequently utilized codons are Leu (UUA), Phe (UUU), and Ile (AUU), as illustrated in Figure 2. It is worth mentioning that, with the exception of Microtendipes bimaculatus (Clade3), the AGG codon (encoding Ser) is absent in the other mitogenomes. Similarly, the CUG codon (encoding Leu) is absent in all mitogenomes except for Microtendipes robustu and Microtendipes wuyiensis. After excluding termination codons, the relative synonymous codon usage (RSCU) was calculated and summarized for the Microtendipes species, as shown in Figure 3. Additionally, we evaluated a suite of codon usage bias metrics, including the Codon Adaptation Index (CAI), Frequency of Optimal Codons (FOP), Effective Number of Codons (ENC), and GC content at the third codon position (GC3), which are detailed in Tables S3–S7. No significant differences were observed among the Microtendipes species in these metrics.

3.4. Phylogeny of Chironominae

In our analysis, all the generic complexes, such as the Polypedilum complex (excluding Microtendipes), Chironomus complex, and Harnischia complex, formed monophyletic groups (Figure 5 and Supplementary Trees S16-S23). The genus Stenochironomus exhibited instability in our phylogenetic trees, being either sister to other Chironominae species or to species within the Polypedilum complex. Similarly, the phylogenetic position of Microtendipes within Chironominae was unstable. Three primary topologies were constructed: (1) Tanytarsus, (Polypedilum complex, (Microtendipes, (Harnischia complex + Chironomus complex))), (2) Tanytarsus, ((Chironomus complex + Harnischia complex), (Polypedilum complex, Microtendipes)), and (3) Tanytarsus, (Microtendipes, (Polypedilum complex, (Chironomus complex + Harnischia complex))), with the last one resembling that proposed by Sæther (1977).

Figure 4. Phylogenetic tree of Chironominae species based on Mitogenome sequences. Values- up on branch present bootstrap supports (%) and Values-down are Bayesian posterior probabilities (PP).

Figure 4. Phylogenetic tree of Chironominae species based on Mitogenome sequences. Values- up on branch present bootstrap supports (%) and Values-down are Bayesian posterior probabilities (PP).

4. Discussion

4.1. 658 bp COI Based Barcode vs Mitogenome-Based Barcode

The 658bp COI barcode has been instrumental in advancing our understanding of species diversity by leveraging its high level of sequence variation and ease of amplification. However, as our understanding of species complexity deepens, the limitations of this barcode have become increasingly apparent. One potential problem exposed by DNA barcodes is the excessive division of species especially in some super diverse genus, e.g. Tanytarsus and Polypedilum [30,46]. Whereas in Song et al. (2023) Microtendipes bimaculatus Song et Qi, with deep intraspecific genetic divergence and forming three paraphyletic clades in phylogenic tree [24]. It is difficult to judge whether there are potential cryptic species in this species group based on the overall situation at that time. While in this study, four specimens mitogenomes from the three clades were used, for Microtendipes bimaculatus-clade 1 (MZ981734); Microtendipes bimaculatus-clade 2 (PP966948); Microtendipes bimaculatus-clade 3 (PP966952 and PP966953), with the genetic distance up to 10.6%. We try to compare the genetic difference between different protein coding genes. The genetic distances were reassessed individually based on 13 protein-coding genes (PCGs), as detailed in Supp. table S8 to S21. The intraspecific distances varied considerably, spanning from 5.3% for NAD4L, 7.0% for NAD1, 7.2% for NAD2, 8.8% for NAD4, 8.2% for NAD5, 9.0% for NAD3, 9.3% for COI, 9.5% for COII, 9.6% for ATP6, 10.0% for COB, 10.0% for NAD6, 10.1% for COIII, to 13.8% for ATP8. Judging from distance-threshold methods, such as insect species proposes a threshold of 2–3%; however, genus-specific thresholds have been suggested for some Chironomidae including 4–5% for Tanytarsus [46] and 5–8% for Polypedilum [30]. Such single gene-based barcodes might mislead when closely related species deep intraspecific splits [46,47]. We also constructed phylogenetic trees to infer the relationships among Microtendipes species separately, using each coding gene, the concatenated 13-protein coding mitogenome genes, and all mitogenome sequences to infer their phylogeny (Supplementary Trees S1-S15). Only the trees based on nad3 and nad4 supported the monophyly of Microtendipes bimaculatus, albeit with low bootstrap values. Therefore, they should be considered part of the Microtendipes bimaculatus species complex, as defined by Song et al. (2023), or sibling species that are difficult to distinguish morphologically [24]. The other divergent species, Microtendipes robustus (PP966950 versus PP966954), formed robust relationships. Further analysis of the variance among different genes was conducted (Figure 4). The pairwise genetic distances of 13 PCGs across nine specimens were computed. For the majority of PCGs, the genetic distances remained statistically indistinguishable. However, two notable exceptions were observed: the genetic distance of ATP8 was significantly higher compared to that of all other PCGs, whereas the genetic distance of nad4l was markedly lower than that of its counterparts.

While the 658bp COI barcode has served as a valuable foundation for species identification and biodiversity research, the mitogenome-based barcode offers a promising alternative that may prove to be more suitable for certain applications. The mitogenome-based barcode’s primary advantage lies in its ability to harness the vast genetic information contained within the mitochondrial genome. This allows for a more nuanced and accurate assessment of species relationships. Furthermore, the mitogenome-based barcode’s comprehensive nature enables researchers to delve deeper into the evolutionary history and population dynamics of species, providing insights that can inform conservation efforts and management strategies.

Figure 4. The genetic distance of PCGs between 9 specimens. CO1_658 represents the 658 bp COI barcode.

Figure 4. The genetic distance of PCGs between 9 specimens. CO1_658 represents the 658 bp COI barcode.

4.2. Phylogenetic Position of Microtendipes

In our analysis, all the generic complexes, such as the Polypedilum complex (excluding Microtendipes), Chironomus complex, and Harnischia complex, formed monophyletic groups (Figure 5 and Supplementary Trees S16-S23). The genus Stenochironomus exhibited instability in our phylogenetic trees, being either sister to other Chironominae species or to species within the Polypedilum complex. Similarly, the phylogenetic position of Microtendipes within Chironominae was unstable. Three primary topologies were constructed: (1) Tanytarsus, (Polypedilum complex, (Microtendipes, (Harnischia complex + Chironomus complex))), (2) Tanytarsus, ((Chironomus complex + Harnischia complex), (Polypedilum complex, Microtendipes)), and (3) Tanytarsus, (Microtendipes, (Polypedilum complex, (Chironomus complex + Harnischia complex))), with the last one resembling that proposed by Sæther (1977) [26]. This similarity suggests that Microtendipes and its related taxa may constitute a generic complex or potentially represent the fourth tribe within this subfamily. Dataset (iii) proposed a relationship of (Chironomus complex (Polypedilum complex (Tanytarsus, Microtendipes))), which does not align ideally with current taxonomy. Specifically, the positioning of Tanytarsus within the Chironomini tribe is incongruent, whereas the similar relationship between Microtendipes and Tanytarsus echoes observations made by Cranston et al. (2012) [28]. Importantly, most results do not support the inclusion of Microtendipes within the Polypedilum complex, thereby refuting the categorization proposed by Sasa (1989) [27]. Given the absence of mitochondrial genome data for Pseudochironomini, we are currently unable to conclusively determine the genus’s taxonomic status. Nonetheless, we are committed to expanding our study by incorporating additional taxa from Tanytarsini and Pseudochironomini, aiming to unravel the precise phylogenetic positioning of Microtendipes with greater accuracy.

4.3. Phylogeny Within Microtendipes

The phylogenetic relationships within the genus Microtendipes have persisted as an enigmatic puzzle, particularly concerning the subdivision into species groups. Two distinct species groups, the pedellus group and the rydalensis group, were established based on larval stages, utilizing two pivotal characters: the number of premandibular teeth and the shape of a median tooth on the mentum [48]. According to this diagnostic criterion, M. baishanzuensis is classified within the rydalensis group, whereas Microtendipes tuberosus and Microtendipes robustus are assigned to the pedellus group (unpublished data). However, it is evident that this classification system does not align seamlessly with the phylogenetic tree, as Microtendipes tuberosus occupies a basal position within the genus’s clade, while other species occupy distinct and separate clades (Figure 5). This discrepancy underscores the need for further investigation to clarify the evolutionary relationships and taxonomic structure within the genus Microtendipes.

5. Conclusions

All eight newly sequenced mitogenomes in this study demonstrated a striking congruity in their structural features and nucleotide compositions when juxtaposed against previously published Chironomidae data. The mitogenome-based species delimitation within the genus Microtendipes offers a more precise and nuanced perspective than COI-based barcoding, particularly those that involve substantial intraspecific genetic distances. The comparative analysis of Microtendipes (Diptera: Chironomidae) mitogenomes has yielded invaluable insights into their phylogenetic relationships, which do not support two species group division based on larvae. This analysis also has definitively refuted the placement of Microtendipes within the Polypedilum complex, instead tentatively establishing it as a distinct species group and tribe, thereby advancing our understanding of this intriguing taxonomic entity.