We used 12 karyologically studied larvae of Ch. degelenus for further DNA extraction. DNA was extracted from the larvae and preserved in 96% ethanol using a Diatom™ DNA Prep 100 kit (Izogen Laboratory Ltd., Moscow, Russia) following the manufacturer’s protocol. DNA extraction was performed on vacuum-dried samples without prior homogenization. Samples were incubated in a lysis buffer at a temperature of 55.5 °C for 18 h. After the extraction, the head capsules were retrieved for dry mounting. The resulting DNA solutions were stored at -18 °C. The amplification of the mitochondrial COI gene was conducted using the MasterMix Х5 kit (Dialat Ltd., Moscow).

To amplify the mitochondrial COI gene's barcoding region, primers LCO1490 (forward) (5′-GGTCAACAAATCATAAAGATATTGG-3′) and HCO2198 (reverse) (5′- TAAACTTCAGGGTGACCAAAAAATCA-3′) (Folmer et al. 1994) were used. PCR was performed in a 25-µL reaction volume. The amplification profile consisted of an initial step of 95 °C for 5 min, followed by 70 cycles of 95 °C for 30 s, 50 °C for 30 s, and 72 °C for 1 min 50 s, and finally a 10-min extension step at 72 °C, a final elongation at 72 °C (10 min), and final storage at 4 °C. The resulting PCR products were purified by precipitation in a 0.15 M CH3COONa solution in 90% ethanol and then rinsed with 70% ethanol. The results were visualized using 1.5% agarose gel electrophoresis with ethidium bromide.

Purified PCR products were sequenced in both directions. DNA sequencing of the COI gene was performed according to Sanger using the BigDye Terminator v3.1 commercial kit (ThermoFisher) and an ABI 3130xl genetic analyser (ThermoFisher) at Syntol JSC (Moscow, Russia). The GenBank accession numbers of the eight Ch. degelenus sequences obtained in this study are OQ873429, OQ873430, OQ873431, OQ873432, OQ873433, OQ873434, OQ873435, and OQ873436.

We used DNA data from both GenBank and BOLD databases (Ratnasingham and Hebert 2007) for the species Ch. balatonicus (JN016826.1), Ch. plumosus (JN016830.1, AB740262.1, BSCHI661-17, BSCHI063-11, GBDP44180-19), Ch. agilis 2 (AF192190.1), Ch. usenicus (JN016806.1), Ch. bonus (MZ014021.1), Ch. entis (GBDPC429-14, KJ085531.1), Ch. borokensis (AB74026.1), Ch. muratensis (AF192194.1), Ch. dorsalis (JN887047.1), Ch. salinarius (KC250756.1, KC250745.1), Ch. annularius (AF192189.1), Ch. bernensis (AF192188.1), Ch. commutatus (AF192187.1), Ch. atrella (USDIQ1762-10, USDIQ1587-10, BBDEE080-10, MHCOL377-07, DRYAS21708-15, SMTPM4256-15), Ch. athalassicus (KM995734.1), Ch. hyperboreus (HQ941615.1), Ch. heteropilicornis (MZ450155.1, MK795770.1, CHMNO266-15), Ch. pilicornis (ARCHR026-11, CNQUF171-14), Ch. staegeri (KR754417.1), Ch. frommeri (KF278235.1), Ch. acutiventris (AF192200.1), Ch. heterodentatus (AF192199.1), Ch. melanescens (AF192204.1), Ch. sokolovae (MW471100.1), Ch. luridus (AF192203.1), Ch. pseudothummi (AF192205.1), Ch. riparius (LN894974.1), Ch. piger (AF192202.1) and Pagastiella orophila (JN265047.1).

The alignment of COI sequences was conducted using MUSCLE with a genetic code for “invertebrate mitochondrial” packaged in MEGA 6 (Tamura et al. 2013). The pairwise sequence distances (Table 2), consisting of the estimated number of base substitutions per site, were calculated using MEGA 6 and the K2P model (Kimura 1980). The analysis involved 51 nucleotide sequences. The codon positions included were the 1st, 2nd, 3rd, and Noncoding. All positions containing gaps and missing data were eliminated. The final data set contained a total of 579 nucleotide positions.

We conducted the estimation of phylogenetic relationships in BEAST V1.10.4 (Suchard et al. 2018) by the Bayesian Markov-chain Monte-Carlo (MCMC) method, using the GTR+G+I substitution model as selected in MEGA 6. The determination of the appropriate model was performed in MEGA 6 (Tamura et al. 2013); the strict clock was used as a clock model and the Yule process as a speciation model. We ran MCMC for 10.000.000 iterations and 1000 iterations of burn-in. Our analysis involved 51 nucleotide sequences, and we eliminated all positions with less than 95% site coverage. There were 579 nucleotide positions in the final data set. We used the COI sequence of P. orophila (Genbank accession number JN265047.1) as an outgroup.

We also tried to get average estimates of divergence time between different branches and clusters that appear on the obtained phylogenetic tree (Figs 3, 4). The age of the most recent common ancestors (TMRCAs) for DNA clades was estimated in BEAST V1.10.4 (Suchard et al. 2018) by the MCMC method, using the GTR+G+I substitution model as selected in MEGA 6. We used a strict clock as a clock model and a constant size as a coalescent model, with the calibration point assumed by Cranston et al. (2012). The time estimate of 36 million years ago (Mya) for the root node of the divergence between P. orophila and all Chironomus species was used as a first calibration point. Also, we use a time estimate of 5.75–3.43 (standard deviations 1.25 and 0.76, respectively) Mya as a second calibration point for the node of the divergence of the species in the Ch. plumosus group of sibling species (Karmokov 2022). We ran MCMC for 10.000.000 iterations and 1000 iterations of burn-in. Tracer v1.7.1 was used to examine the BEAST log file and ESSs for each parameter, which were all > 200.

According to recent studies (Jamnongluk et al. 2003; Stevens et al. 2006), the COI gene sequence diverges at rates ranging from 1.5% to 2.3% per 1 Mya in insects. Under the proposed partitioning scheme and substitution model chosen using Bayes factors, Papadopoulou et al. (2010) found a divergence rate for the COI gene in the study of tenebrionid beetles of 3.54% per 1 Mya (2.69% when paired with the 16S rRNA gene). In this investigation, we used the three frequently accepted mutation rates of 1.5%, 2.3%, and 3.54% to calculate the divergence time. We determined TMRCAs for the nodes 1 through to 5 of the phylogenetic tree (Fig. 3). The results are presented in Table 3.

The diploid number of chromosomes in the Ch. degelenus karyotype is 2n = 8, and the chromosome arm combinations are AB, CD, EF, and G (the “thummi” cytocomplex) (Fig. 3). Chromosome AB is metacentric, CD and EF are submetacentric, and G is telocentric. The centromeric bands are prominent. There are two nucleoli: one of them is located on the centromeric end of the arm G, and another is in the middle of the arm B. Two Balbiani Rings are developed on the arm G. In arm G, the degree of development of the Balbiani ring and the degree of homologue conjugation in different larvae varied.

Previously, Seisebaev et al. (2001) described 7 banding sequences in the Ch. degelenus banding sequence pool. In our study, all of those sequences are present, and one new banding sequence, degC2.2, has been found for the first time, providing eight banding sequences in the Caucasian population.

We should note that after we compared the obtained Ch. degelenus COI sequences with those in the databases (using the BLAST (The Basic Local Alignment Search Tool) algorithm in the NSBI database), we found that the closest species to Ch. degelenus is Ch. atrella (Townes, 1945). Looking through the available data in the literature, we see that the karyotypes of these species are surprisingly similar. The patterns of the banding sequences of Ch. degelenus have almost identical analogues among the banding sequences of Ch. atrella. For our comparison we used the photomaps of Ch. atrella from Martin et al. (2006). In our opinion, the banding sequences in all chromosome arms, except arm G, are highly similar if not identical. However, the mapping in most parts of the arms is different. To illustrate this, we provide photomaps of Ch. atrella with the mapping of Martin et al. (2006) for each arm, along with photomaps of Ch. degelenus with the mapping of Seisebaev et al. (2001) (Fig. 4). The mapping of the original samples by Seisebaev et al. (2001) was followed in the mapping of Caucasian samples.

Arm A has one banding sequence degA1 (Figs 3, 4). The banding patterns of sequence degA1.1 are completely identical to atrA1.1 except for one small feature in the vicinity of the centromere. In other species of Chironomus, the centromeric band of this arm is normally represented by a thick band (K), located immediately after the thick-band group 19ef (Keyl, 1962). However, in Ch. atrella, there is a significant region containing at least three pale bands between 19ef and K (Martin et al. 2006). The reason for such a structural modification is not clear. Such a region with pale bands is also present in Ch. degelenus samples from Kazakhstan but absent in Caucasian ones.

Arm B was monomorphic with banding sequence degB1.1 (Fig. 5). The banding patterns of the sequence is highly similar to banding sequence in arm B that Martin et al. (2001) signed as atrB2+B3.

Arm С has two banding sequences: degC1 and degC2 (Fig. 6). The banding sequence degC2 is new for the species and described for the first time. The banding sequence degC2 (occurrence: 0.708) and genotypic combination degC2.2 (occurrence: 0.542) were predominant in the Eastern Caucasus population o. The sequence degC1 (occurrence: 0.292) and genotypic combination degC1.1 (occurrence: 0.125) were found much less often in the population from the Eastern Caucasus. The heterozygote degC1.2 was also observed with a low frequency (occurrence: 0.333). The sequence degC2 differs from degC1 by simple inversion in the area 8a-11c 4f-a.

The banding patterns of the sequence degC1 (Fig. 6D) is identical to the sequence atrC2 (Fig. 6C) and the sequence degC2 (Fig. 6A) is identical to the sequence atrC1 (Fig. 6B).

Arm D also was monomorphic with banding sequence degD1.1 (Fig. 7). The banding patterns of the sequence degD1 (Fig. 7A, C) is identical to the sequence atrD2 (Fig. 7B).

Arm E was monomorphic with banding sequence degE1.1 (Fig. 8). The banding patterns of the sequence degE1 (Fig. 8A, C) is identical to the sequence atrE1 (Fig. 8B). It is the only arm where the chromosome mapping of Martin et al. (2006) for Ch. atrella and Seisebaev et al. (2001) for Ch. degelenus correspond.

Arm F was monomorphic with a banding sequence degF1.1 (Fig. 9). The banding patterns of the sequence degF1 (Fig. 9A, B, D) are identical to those of the sequence atrF1 (Fig. 9C). Seisebaev et al. (2001) did not note the presence of a nucleolus or Balbiani rings in the arm, but in Caucasian larvae, one can clearly see that in the arm there is some structure resembling a nucleolus. Interestingly, Ch. atrella also has a nucleolus in the same region (Martin et al. 2006). The sequence degF1 is still not mapped, while the sequence atrF1 is mapped by Martin et al. (2006), and probably one can use this mapping for sequence degF1.

Arm G was monomorphic with banding sequence degG1.1 (Fig. 3, 10). This is the only arm in the karyotype where chromosome structure is clearly different in both species.

The level of inversion polymorphism in the Caucasian Ch. degelenus population is quite low. Chromosomal polymorphism was observed only in arm C. The percentage of heterozygous larvae in the population of the Eastern Caucasus is low (33.3%). The average number of heterozygous inversions per larva is also relatively low (0.333).

Overall, we successfully obtained 8 sequences of Ch. degelenus from 12 Caucasian larvae. We found six haplotypes among the obtained sequences. Two haplotypes occur twice each among the obtained sequences: the first haplotype: OQ873433, OQ873431, and the second haplotype: OQ873436, OQ873429, respectively. The CG content and base composition of the obtained sequences are provided in Table 1. This is the first DNA data obtained for Ch. degelenus.

As we found that the closest related species to Ch. degelenus was Ch. atrella, we decided to look closer at the Ch. atrella DNA data available in the BOLD database. In general, this species is predominantly found in North America. There are a total of 559 COI gene sequences for Ch. atrella in the database. They formed two BINs (clusters) with specimens from two countries: Canada (548 sequences) and the USA (11 sequences). The Barcode Index Number System (BINs) in BOLD is an online framework that clusters barcode sequences algorithmically, generating a web page for each cluster. Since clusters show high concordance with species, this system can be used to verify species identifications as well as document diversity when taxonomic information is lacking.

In our analysis, we decided to use the sequences from both BINs and take three haplotypes from each one. The BIN with the largest number of sequences was assigned as the 1st BIN, and the other one as the 2nd BIN. We should note that almost all of the COI sequences available in the BOLD database of Ch. atrella were identified by the BOLD Identification System (BOLD-IDS), i.e., through COI gene barcodes without karyological analysis.

The resulting phylogenetic tree (Fig. 11) shows a very complex pattern, each having relatively high probabilities. Some of these clusters correspond well with the groups of closely related species, such as the Ch. plumosus, Ch. pilicornis, Ch. piger, and so on. The most interesting clusters for us are those that contain sequences of Ch. degelenus and Ch. atrella.

We conditionally named them the Ch. degelenus cluster, the Ch. atrella 1st BIN cluster, and the Ch. atrella 2nd BIN cluster.

The Ch. degelenus cluster (Fig. 11) is formed by Ch. degelenus sequences from Caucasian specimens. The closest cluster to this is the Ch. atrella 1st BIN cluster, which is formed by sequences of Ch. atrella from 1st BIN. Next, the closest to these two clusters is a sequence of Ch. athalassicus Cannings, 1975. There are a total of 523 COI gene sequences for Ch. athalassicus in the BOLD database. They formed one BIN (cluster) with specimens from two countries: Canada (522 sequences) and the USA (1 sequence).

The sequences of Ch. atrella from the 2nd BIN form their own separate cluster, which at some distance is joined by Ch. hyperboreus Staeger, 1845. There are a total of 429 COI gene sequences for Ch. hyperboreus in the BOLD database. They formed two BINs (clusters) with specimens from three countries: Canada (317 sequences), Greenland (109 sequences), and Norway (3 sequences).

Calculated pairwise sequence COI gene distances (Fig. 12, Table 2,), consisting of the estimated number of base substitutions per site using the K2P model (Kimura 1980), reveal an interesting pattern. Proulx et al. (2013) demonstrated that intraspecific K2P distances for Chironomus species, based on the COI gene ranged from 0% to 3%. They did this by using genetic, morphological, and karyological data to distinguish Chironomus species from Canada. These values could potentially be used as a guide in the current study to discriminate between different Chironomus species. However, information from the COI gene must be supplemented by other methods.

In the study, the distances between the sequences within the Ch. degelenus cluster, the Ch. atrella 1st BIN cluster, and the Ch. atrella 2nd BIN cluster (Table 2) are below the 3% threshold. The distances within the Ch. degelenus cluster are very small, ranging from 0% to 1.65% (average: 0.9%). The distances within the Ch. atrella 1st BIN cluster are even smaller, ranging from 0% to 0.36% (average: 0.2%), while those within the Ch. atrella 2nd BIN cluster vary from 0% to 0.55% (average: 0.5%).

In contrast, the average distances between the various clusters exceed the 3% threshold. The distance between the Ch. degelenus cluster and the Ch. atrella 1st BIN cluster is 5.1%, while the distance between the Ch. degelenus cluster and the Ch. atrella 2nd BIN cluster is 6.6%. Finally, the distance between two Ch. atrella clusters is 6.2%.

The earliest split between the Ch. degelenus cluster and Ch. atrella clusters, according to our calculations, occurred between 3.72-2.67 Mya (Fig. 11; Table 3, node 3; substitution rates for all earliest and latest estimates in this chapter are 1.5% and 3.54%, respectively). These estimates fall within the Piacenzian epoch (3.60-2.58 Mya), which is the most recent stage of the Pliocene era (5.33-2.58 Mya). The most recent common ancestor of the Ch. degelenus cluster and the Ch. atrella 1st BIN cluster lived 2.05 and 1.47 Mya (Fig. 11; Table 3, node 4). This division took place during the Calabrian epoch of the Pleistocene. The most recent common ancestor of Ch. degelenus cluster lived between 0.504 and 0.356 Mya (Fig. 11; Table 3, node 5), corresponding to the Chibanian stage of the Pleistocene (0.770-0.126 Mya).

In conclusion, the present study has provided valuable insights into the phylogenetic relationships and chromosomal characteristics of Chironomus species. The discovery of Ch. degelenus in the Eastern Caucasus and Russia for the first time is a notable finding, and the low level of polymorphism observed in the Caucasian population of the species underscores the need to explore further genetic diversity of the species. Additionally, the identification of a new banding sequence, degC2, in the banding sequence pool of Ch. degelenus, adds further evidence of chromosomal diversity in this species.

The observed similarity between the polytene chromosome banding patterns of Ch. degelenus and Ch. atrella suggests that these two species may be close relatives. However, the analysis of K2P distances between the clusters of Ch. degelenus and Ch. atrella supports their classification as separate species. The observed picture with the species Ch. degelenus and Ch. atrella could be seen as one of the numerous manifestations of a phenomenon known as “dark taxa”, which can greatly complicate the work of specialists, especially experts in species of the genus Chironomus and closely related genera.

We propose that the emergence of grasslands in the Pliocene created new habitats for Chironomus species, and this, combined with the role of wind in dispersal, may have been a key factor in the divergence of Chironomus species, including Ch. degelenus, Ch. atrella (both BINs), Ch. athalassicus, and Ch. hyperboreus. Further research is needed to understand the mechanisms driving speciation in these species and to evaluate the factors responsible for the observed chromosomal divergence.

Overall, this study provides important new insights into the evolutionary history and genetic diversity of Chironomus species, shedding light on the complex processes that underlie their divergence and speciation.