Molecular Phylogeny of Holarctic Aeshnidae with a Focus on the West Palaearctic and Some Remarks on Its Genera Worldwide (Aeshnidae, Odonata) †

: Aeshnidae Rambur, 1842 are impressive large insects distributed worldwide. Currently, over 500 species are recognized. Nevertheless, the phylogeny of this family is not completely understood. We applied molecular phylogenetic analysis using two popular phylogenetic markers, the mitochondrial COI gene fragment (barcoding sequence) and the nucleic ITS region, containing the ITS1, 5.8S rRNA, and ITS2 sequences. We used available and credible published sequences and 96 newly sequenced specimens. Our analysis involved all West Palaearctic species, all but one genera of the Holarctic Aeshnidae, and most genera worldwide, and is by far the largest molecular study of this family. The topology of all trees created with different algorithms and genes is in favour of the current taxonomic concept, with some remarkable outcomes. Aeshna Fabricius, 1775, was found to be diverged into several branches, especially with respect to the COI gene. Although it appeared not monophyletic in phylogenetic reconstructions based on the ITS region, the analysis of COI and joint analysis suggest its monophyly in the current taxonomical sense, with one notable exception. Aeshna isoceles (Müller, 1767) has fallen out of Aeshna in all analyses, so a new monophyletic genus, Isoaeschna gen. nov. is introduced for it. The genus Brachytron Evans, 1845 tightly clustered with Aeschnophlebia Selys, 1883, Epiaeschna Hagen in Selys, 1883, and Nasiaeschna Selys in Förster, 1900. Thus, we suggest subsuming these four genera under the priority name Brachytron . Tetracanthagyna Selys, 1883 clusters as expected with Brachytron in the ITS tree, but is an independent ancient clade of its own in all COI trees. The genus Polycanthagyna Fraser, 1933 syn. nov. is synonymised to Indaeschna Fraser, 1926. On the species level, we suggest that the American Aeshna septentrionalis Burmeister, 1839 be treated as a subspecies of A. caerulea (Ström, 1783), Aeshna caerulea septentrionalis . We synonymize Gynacantha hyalina Selys, 1882 with Gynacantha subinterrupta Rambur, 1842. Our analysis provides new insights on the tight relationships of the circumboreal species Aeshna juncea and A. subarctica and the intraspecies phylogeny of Aeshna juncea .


Introduction
Aeshnidae Rambur, 1842 is a diverse family of Ansioptera (Figure 1) embracing large and robust dragonflies, which are strong fliers; some are even able to cross oceans from one continent to another.The family is distributed worldwide, and the number of species included is steadily increasing and is currently well above 500 [1].The representatives of Aeshnidae are currently sorted in 54 accepted genera [1].Less than half of the genera and about one-fifth of species is found in the Holarctic.Several classifications for Aeshnidae were proposed mainly on arbitrarily chosen morphological characteristics and wing venation [2][3][4].More selected characters from external morphology have been used in a cladistic analysis of all genera of Aeshnidae known up to 2001 by Natalia von Ellenrieder [5].Studies using molecular methods mostly addressed other families or the entire order of Odonata, the latter including only a small number of Aeshnidae [6][7][8][9][10][11][12].A recent study on the genus Anax Leach in Brewster, 1815 brought some light into the relationship of migration and phylogeny [13].These studies left the phylogenetic relationships between the Holarctic Aeshnidae unresolved.
A molecular genetic study of Aeshnidae that includes a considerable number of specimens is still missing.Therefore, we made such an attempt, focused on two popular markers, the mitochondrial COI gene fragment and the nuclear ITS regioninvolving both published sequences and 96 newly sequenced specimens.We used several algorithms to obtain the most appropriate results.Our analysis comprises nearly all West Palaearctic species and nearly all genera of the Holarctic Aeshnidae and is by far the largest molecular study of this family so far.

Materials
Our molecular analysis included a total of 291 Aeshnidae specimens (Tables 1 and S1).Of them, 96 specimens were sequenced in the course of this study, 27 specimens were sequenced and published by us previously [14], and the sequences of 168 specimens were taken from GenBank [15].The latter were earlier published in [10,16,17], and some other papers (Table 1; for more information, including the PCR numbers mentioned in the figures with trees, consult Table S1).Our analysis involved 28 of 54 currently recognised genera of Aeshnidae, and 17 of 18 genera occurring in the Holarctic region (with Gomphaeschna Selys, 1871 being the only genus not involved).Eleven genera from beyond the Holarctic region were also included, as they were supposed to have relatives in the Holarctic region.For the ITS region, which contains the ITS1, 5.8S rRNA, and ITS2 sequences, the sequences are between 648 and 1006 bp long, depending on the species (Figure 2).For the barcoding fragment of the COI gene, we used two alignments: one with sequences with a length of 632 bp (Figure 3) and the other with shorter sequences, also available in GenBank [15], so that all sequences were trimmed to the same length of the shortest, 341 bp (see Section 3.2.2).We used Orthetrum melania (Selys, 1883) and Neopetalia punctata (Hagen in Selys, 1854) as outgroups for the ITS and COI analyses, respectively.

DNA Extraction and Sequencing
Per specimen, a 1 mm section of a leg was transferred to a tube with 20 µL 0.05 N NaOH and 2 µL 5% Tween 20.This was heated for 15 min at 95 • C and cooled on ice.A volume of 100 µL sterile water was added to the tube and mixed.The amount of 1 to 5 µL of this solution was used in a PCR reaction.For details of PCR and sequencing see [33,34]; for the COI barcoding fragment we used the primers CO1490F (50-GGT CAA ATC ATA AAG ATA TTG G-30) and CO2198R (50-TAA ACT TCA GGG TGA CCA AAA AAT CA-30) and for the ITS region we used the primers Vrain2F (50-CTT TGT ACA CAC CGC CCG TCG CT-30) and 28R1 (50-TGA TAT GCT TAA NTT CAG CGG GT-30).

Reducing Artefacts
All sequences used were checked for plausibility of determination.Each sequence was blasted to check for sequencing contamination and in case of doubts, it was checked to determine that it was not a pseudogene by comparing it with other sequences and translating it to protein.This was also done for some GenBank sequences that appeared in unexpected branches in the tree.Some sequences were revealed to be pseudogenes and were discarded.The alignments with a COI sequence of Austrogynacantha heterogena Tillyard, 1908 from GenBank revealed a strange, close relationship with Aeshna isosceles; therefore, we isolated DNA from an additional specimen of this species we received from the Australian Museum and approved the correctness of the GenBank sequence.
We also eliminated positions in the alignment that were saturated with multiple substitutions with the program Gblocks (Figure S2) [35] and constructed trees with the so-reduced alignment (see below).

Phylogenetic Analysis
Alignments were made with the online version of MAFFT [36] with default settings.The model of DNA evolution that best fit the data was determined with JMODELTEST version 2.1.10[37].Based on the Bayesian information criteria (BIC), the best model was chosen (nst = 6, rates = gamma for both COI and ITS analysis).With this model of evolution, trees were constructed using MRBAYES 3.2.7a[38].The settings were as follows: 10 million generations, a sample frequency of 1000, and a burnin value of 5000 trees.For more detail, see [33,34].
Since the COI-based trees were less resolved than the ITS-based tree, we also tried to eliminate positions in the alignment that were saturated with multiple substitutions with the program Gblocks (Figure S1) [35] and constructed trees with the so-reduced alignment that was 315 bp long.
In addition, both sequences altogether were analysed with StarBeast3 v1.1.7 [39], which is a multi-individual multi-locus species tree estimation program, using Bayesian coalescent analysis, as implemented in the BEAST v2.7.3 package [40].This approach takes into account that sequences do not evolve alone but are always present in some species which may originate from each other by divergence.Xml input files were created in BEAUTI v2.7.4,using the HKY + Γ + I model for both markers.The following settings were used for all analysis: base frequencies: 'empirical'; clock model: 'Strict clock Clock.rate:1'; TreePrior: 'Yule Model'; popMean: Log Normal with M: −5 ans S: 1.2; clockRates: 'Exponential'.The analyses were run on BEAST software.Analyses were run for 10 million generations, sampling every 5000th generation.Tracer v. 1.7.1 [41] was used for examining the effective sample size (ESS) for parameters and determining the burnin.Trees and posterior probabilities were summarized using TreeAnnotator v. 2.7.3 and shown on the maximum clade credibility tree with median heights, with a posterior probability limit = 0.5 and burnin percentage = 10.The trees were drawn in FigTree v.1.4.4 (http: //tree.bio.ed.ac.uk/software/figtree/, accessed on 1 June 2023).

Haplotype Network Analysis
Haplotype networks were built based on the COI alignment using POPART 1.7 [42] software with the TCS network inference method.A haplotype network is the sum of the shortest evolutionary pathways between the current haplotypes via subsequent mutations that connect the current DNA molecules via putative intermediate molecules.

Analysis Based on the ITS Region
We begin our analysis with the more conserved gene fragment of the nuclear ribosomal RNA region including the intervening ITS region: ITS1, the 5.8S rRNA coding sequence, and ITS2.A phylogenetic tree reconstructed on the base of the ITS regions by a Bayesian approach is provided in Figure 2. The left part of this figure presents the tree as a phylogram, in which lengths of branches are proportional to the number of accumulated nucleotide substitutions, which depends on both the time of divergence and the rate of molecular evolution.The right part of Figure 2 is the same mirrored tree as a cladogram, which shows only the tree topology.The cladogram is added to visualise divergence and the clustering of sequences differing in few substitutions only, which is not seen in the phylogram (on the left) because of too-short branches.
In general, the tree based on the ITS region corresponds well to the current taxonomy of Aeshnidae and shows most of the current genera as monophyletic, with one notable exception of the largest and most familiar genus, Aeshna Fabricius, 1775.In this tree, 16 clades can be recognized.
Clade 1 (juncea-clade) comprises two species, Aeshna juncea (Linnaeus, 1758) and Aeshna subarctica Walker, 1908.Surprisingly, A. juncea was represented by two clusters (not well seen in the phylogram on the left side of Figure 2 but visualised in the cladogram of its right side), one with specimens from the Caucasus and Transcaucasia region and the other including the rest of specimens from elsewhere.Unfortunately, no American representative of A. juncea was available for the ITS region analysis.A. juncea and A. subarctica appeared as sister species, as expected.
Clade 2 (grandis-clade) is composed of Aeshna grandis (Linnaeus, 1758), which is the type species of the genus Aeshna Cowely (1934), Aeshna viridis Eversmann, 1836, Aeshna serrata Hagen, 1856, and Aeshna crenata Hagen, 1856.Strikingly, this clade of these four species diverges in only two monophyletic clusters, that for A. crenata and that for, altogether, A. serrata, A. grandis and A. viridis (Figure 2, see the cladogram on the right).Curiously, the three latter species, so different in appearance and ecology, look in the tree as if they were a single species, while A. crenata, representing the other cluster, is a species strongly resembling A. serrata in appearance and even sometimes confused with it.2, see the cladogram on the right side).This dichotomy has been recognized before [14,44], and is paralleled by the above-discussed divergence of A. juncea from the Caucasus and Transcaucasia versus from elsewhere.The American members of the 'umbrosa group' ('paddle-tails') were not available for the ITS analysis.
Clade 4 includes Aeshna mixta Latreille, 1805 and Aeshna soneharai (Asahina, 1988) as very closely related sister taxa.Actually, their ITS region differs in just two substitutions in the ITS region.It is noteworthy that Onishko et al. [23], who raised the taxon A. soneharai to species level based on external characters, behaviour, sympatric occurrence, and differences in the mitochondrial COII sequences, also sequenced the ITS2 spacer and found it to be identical to A. mixta.We now have sequenced the broader ITS region and found that Aeshna soneharai differs in three nucleotide substitutions from A. mixta.
Clade 9 is distinctly separated from all the above-mentioned clades and includes only one species: Aeshna isoceles (Müller, 1767), with all its sequences identical, although our analysis contains specimens comprising the whole geographical distribution of this species, including Europe, the Near East, West Asia, and North Africa.
It is noteworthy that the cluster uniting clade 9 with clades 1-8 has the negligible support of a 0.56 posterior probability value, which provides a strong argument against inclusion of A. isoceles into the genus Aeshna.The cluster uniting clades 1-8 has the highest support, 1.0, and could be considered as one genus named Aeshna, but this, at the same time, would imply synonymization of the genera Anaciaeschna, Polycanthagyna, Andaeschna, and Pinheyschna Peters et Theischinger, 2011 with Aeshna.This view is, however, not supported by the COI analysis (see below).
Clade 10 represents the genus Anax.The clade includes Anax ephippiger (Burmeister, 1839) and makes an additional genus as Hemianax Selys, 1883 unnecessary, which is in line with previous studies, for example the most recent one by [13].Anax imperator Leach in Brewster, 1815 and Anax parthenope (Selys, 1839) are well separated, while A. parthenope (unfortunately represented in the ITS tree by a single specimen only) and Anax julius Brauer, 1865 are not separated.
Clade 11 is represented in our analysis only by one specimen and species, Rhionaeschna diffinis (Rambur, 1842) from Chile.Few members of this genus reach the Holarctic in southern North America.A better relation of this genus to the others can be seen in the StarBeast analysis of the combined genes (see below).
Clade 12 is composed of three genera: Gynacantha Rambur, 1842, Austrogynacantha Tillyard, 1908 and Heliaeschna Selys, 1882.Members of the last two genera are not present in the Holarctic.
Clades 13, 14 and 15 are an interesting complex of genera.All these genera have very strong supports on the tree and could be made clades of their own, so our subdivision of this complex into clades is rather arbitrary.Clade [5], as containing the first two of its three groups.These three clades may be assumed as the subfamily Brachytroninae, as suggested previously [4].Besides the morphological similarity, the members of this subfamily share also similarities in behaviour and habitat selection, with most of them preferring shady stream sections or marshes.
Clade 16 is represented by Sarasaeschna McLachlan, 1896 only (the node uniting it with clade 15 is scarcely supported, the posterior probability being 0.73).

Analysis Based on the COI Gene
For the analysis based on the COI gene fragment, we attempted several options.First, we reconstructed phylogenies based on long (632 bp) and short (339 bp) fragments of the COI gene.For the purpose of this study, we sequenced the long fragment.The short fragment was naturally less informative but had an advantage of having much better representation in GenBank [15], so that we were able to include many more species.We also attempted a Gblocks analysis, which removes positions in the alignment that are saturated by substitutions or poorly aligned.These three analyses yielded similar results, and although the long fragment tree included fewer species, it was highly representative for the Holarctic because of our efforts to de novo sequence relevant specimens.Therefore, we will describe below the phylogenetic tree reconstructed on the base of the long fragment (Figure 3).

Analysis Based on Long COI Fragment
Compared to the ITS tree, the COI tree (Figure 3) contained much more basic clades, including smaller number of sequences.Furthermore, the COI sequence did not resolve the phylogenetic relationships of quite a number of those basic clades, thus revealing polytomy.Nevertheless, grouping of genera in this tree appeared very interesting.The uppermost main node of the tree (Figure 3), with a weak support of 0.68, corresponds to the ITS clade 3 and can be called the cyanea-clade.It contains three well supported branches.The first of them is composed of A. cyanea and A. vercanica; the second includes the North American members of the 'umbrosa group' ('paddle-tails'), Aeshna umbrosa Walker, 1908, Aeshna constricta Say, 1840, and Aeshna palmata Hagen, 1856, the sequences of which were available in GenBank (2023).The third branch includes the Asian A. petalura.
Specimens of A. constricta Say, 1840 and A. palmata Hagen, 1856 are interspersed in the tree (Figure 3), as if they were the same species.The haplotype network (Figure 4) shows that they share the most common allele of the studied COI fragment.and A. petal.Like in the ITS analysis, specimens of A. cyanea from North Africa and Europe cluster together.However, specimens of A. cyanea from the Caucasus and Transcaucasia region on one hand and from the rest of the range on the other had no longer form sister clades as in the ITS tree (Figure 2).Instead, the Caucasian plus Transcaucasian specimens radiate from the base of the A. cyanea cluster, while the European plus African cluster is now internal, as a sprouting among them.This can be interpreted as the species A. cyanea having originated in the Caucasian area and then one of its lineages having spread to the west and occupied vast European and North African territories.A. vercanica and A. cyanea are sister branches in the COI tree (Figure 3), but the node of A. cyanea has a weak support of 0.55.
The next clade, with the highest possible support of 1.0 (the grandis-clade) corresponds to the ITS clade 2.Besides the West Palaearctic A. grandis, A. viridis, A. serrata, and A. crenata, it also includes the North American Aeshna interrupta Walker, 1908, Aeshna eremita Scudder, 1866, Aeshna canadensis Walker, 1908, Aeshna verticalis Hagen, 1861, and Aeshna tuberculifera Walker, 1908 in the COI tree (Figure 3).The species A. grandis and A. viridis form clusters of their own, having rather recently diverged but having the maximum support of 1.0, and are not united with A. serrata as they were in the ITS tree (Figure 2).
Two putative cases of introgression between A. serrata and A. crenata were detected.Two specimens, from Finland and the Ural Mountains (Russia), were identified as A. crenata but had a COI sequence identical to A. serrata.Since their identification by morphological means was unequivocal (the specimen from Ural was collected and examined by one of us), we may suggest that this was an old mitochondrial introgression from A. serrata into A. crenata rather than recent hybridization.
The next well supported clade includes the Eurasian species A. caerulea and the North American species Aeshna septentrionalis Burmeister, 1839 and A. sitchensis Hagen, 1861.The node uniting this clade to the previous one is not supported (0.73), so should not be taken into account.The branch of A. septentrionalis is not a sister one to any other species, but appeared as an inner branch inside A. caerulea.This is also well illustrated by the haplotype network constructed for this clade (Figure 5), where alleles revealed in the two Canadian species appeared to independently originate from that found in a specimen from Austria.So, our data rule out the species level of A. septentrionalis.Therefore, we synonymize it with A. caerulea at the species level as suggested before [46][47][48], downgrading it to the subspecies of the latter.The next large clade (juncea-clade) in the COI analysis corresponds to the ITS clade 1 and comprises the same two species, A. juncea and A. subarctica.Like in the ITS analysis, there are two clusters, but their content is striking.One of them comprises specimens of A. juncea originating from the West Palaearctic, up to West Siberia to the east.The specimens from the Caucasus region (Armenia, Georgia, North Caucasus) are no longer separated but are interspersed with European specimens.The second cluster contains specimens of A. juncea originating from the eastern half of Eurasia (Pakistan and the Russian Far East), and from Canada (North America).This result is in agreement with that reported by Kohli et al. [11], who found common COI haplotypes of North American and East Asian (Japan and China) specimens of A. juncea, which were different from those of European specimens.Unexpectedly, the second cluster contains, also, all specimens of A. subarctica.Thus, the structure of the juncea-clade in the COI tree contradicts not only that in the ITS tree but also the long-established morphological systematics.Obviously, we faced a case wherein mitochondria exhibit a phylogeny of their own, discordant to that of nuclear sequences and that resulting from morphological data.Such cases frequently appear in Odonata [20,49,50].This also underlines the close relationship between A. juncea and A. subarctica.The COI haplotype network of the two species shows, from a different aspect, the same pattern, in which the Far Eastern and North American specimens of A. juncea are more separated from the West Palaearctic A. juncea than from A. subarctica (Figure 6).The next clade in the COI tree corresponds to the ITS clade 4 but includes A. affinis as its first divergence, as expected.In the COI tree, A. mixta and A. soneharai are not two sister branches, as in the ITS tree, but A. soneharai appeared as an inner branch inside A. mixta.At the same time, in the haplotype network (Figure 7), these species are independent branches.The Kimura 2-parameter distance (for details see [33]) between A. mixta and A. soneharai is small (0.02), suggesting very recent separation of these taxa (Figure 7).For sequences with so few substitutions, a haplotype network is a more adequate representation than a phylogenetic tree, since correct phylogenetic analysis demands a considerable signal from many substitutions.The next large cluster in the COI tree is not supported.Inside this, a well supported branch represents the genus Boyeria, with four Holarctic species; the next one contains two species of Planaeschna, P. cucphuongensis Karube, 1999 and P. risi Asahina, 1981.The third species of this genus, P. ishigakiana Asahina, 1951, fell aside in a polytomic cluster, with a negligible support of 0.62, with the two previous species and Boyeria.The next well supported branch is represented by Polycanthagyna erythromelas (McLachlan, 1896) and Indaeschna grubaueri (Förster, 1904), tightly clustering together with the highest possible support, 1.00.The same clustering of these two species from different genera is supported also by the combined anaylsis by StarBeast (see below).The following group of eight sequences represents the genera Caliaeschna, Sarasaeschna, Andaeschna, and Periaeschna but it does not form a clade, as neither of its nodes is supported.Even representatives of the same species do not show significant clustering.For instance, all the four specimens of Caliaeschna microstigma (Schneider, 1845) cluster to each other with supports not higher than 0.76, with a specimen from Azerbaijan distant from the other three (including the second specimen from Azerbaijan).The two specimens of Periaeschna magdalena Martin, 1909 do not cluster with each other at all.The genera Cephalaeschna, Periaeschna, and Planaeschna, which clustered with Caliaeschna in the ITS tree, are now sorted apart.The haplotype analysis also showed that Caliaeschna is not closely related to Cephalaeschna, Planaeschna, Sarasaeschna, and Remartinia Navás, 1911 (Figure 8).Thus, Caliaeschna cannot be united with Cephalaeschna, as discussed previously [51].Our specimens of Caliaschna microstigma (the only species of its genus) originate from all over its range, from Montenegro to Azerbaijan.Their cluster is not supported, although the sequences do not differ much, as suggested by the haplotypic network (Figure 8).In particular, a specimen from Balakan District of Azerbaijan (sequence 16161) represents a local population characterised by strongly reduced antehumeral stripes, not paralleled by shrinkage of other pale markings [52].At the same time, the specimen from Ordubad District of Azerbaijan (sequence 16160) has normal, not reduced antehumeral stripes.The reduced stripe was observed in the neighbouring Georgia [53], eastern parts of Turkey (unpublished), while in some populations of Dagestan (Russia) and SE Turkey, individuals with both stripe versions fly together (unpublished).Both our analyses, of ITS (Figure 2) and COI (Figure 3), did not reveal any divergence of the Balakan specimen from those from elsewhere; its COI sequence is most close to the specimen from Montenegro.This fact suggests that the reduction of the antehumeral stripe does not manifest a special Caucasian taxon, even of a subspecific rank.
The next three isolated branches are composed of (i) three specimens of Basiaeschna janata (Say, 1840), representing the monotypical American genus Basiaeschna Selys, 1883, (ii) Remartinia liteipennis (Burmeister, 1839) and Coryphaeschna adnexa (Hagen, 1861), but the support of this branch of 0.57 is negligible; (iii) two specimens of Cephalaeschna risi Asahina, 1851, the only representative of its large genus in our analysis.
The next large and well supported clade includes only the genus Gynacantha, represented in our analysis by quite a number of species.Two major, well supported branches can be recognized in this clade, one representing the African members, the other the Asian members of the genus (as seen from country annotations at the species).A corresponding topology is seen in the haplotype network, with the species Gynacantha bispina Rambur, 1842 from Mauritius placed in between (Figure 9).The taxonomy of this genus in Africa is problematic [54].Looking deeper on the species level, a high similarity of Gynacantha congolica and G. manderica was seen, and the difference in the haplotype tree was beyond the species level (Figure 9); however, we did not have the corresponding ITS sequences to definitively synonymize them.At the same time we propose to synonymize the Asian Gynacantha hyalina Selys, 1882 with Gynacantha subinterrupta Rambur, 1842; this is supported by the ITS (Figure 2) and COI (Figure 3) trees and the haplotype network (Figure 9).The next clade includes A. isoceles only, so corresponds to the ITS clade 9. We have analysed specimens from distant geographical regions comprising the whole range of the species and found it highly homogeneous.There is no place for subspecies like A. isoceles antehumeralis Schmidt, 1950 (alternatively Anaciaeschna isoceles antehumerlis) as suggested by Schmidt [43].Both markers analysed did not reveal any relationship of this species to Anaciaeschna, as supposed previously [43,55].The relationship to Austrogynacantha (Figure 3), as also previously illustrated by Carle [9]'s supplement (COI tree), seems to be an artefact, since the node is weakly supported (the posterior probability being 0.73).It is noteworthy that this node is no longer present in the Gblock tree (Figure S3).Also, the haplotype network shows a far separation of both genera (Figure 10).Furthermore, this relation was not indicated by the ITS analysis.The clustering of the branch formed by a single specimen of Staurophlebia reticulata (Burmeister, 1839) with A. isoceles and Austrogynacantha is insignificant (0.73).While the ITS tree showed well the subfamily Brachytroninae, represented by three clades 13-15 (see above), in the COI tree it is no longer traced, but split in as many as six branches (Figure 3).
A robust clade with the maximum support is composed of the genera Brachytron, Aeschnophlebia, and Epiaeschna Hagen in Selys, 1883.The West Palaearctic genus Brachytron consists of a single species, Brachytron pratense (Müller, 1764), which is represented in our analysis by specimens from throughout its geographical range.The haplotype network (Figure 11) also suggests a very close proximity of representatives of the above genera.The genus Tetracanthagyna, which clustered in the ITS analysis with the genus Brachytron, is now far outside this cluster, forming the most basic clade of the Aeshnidae tree (Figure 3).Such different sorting of the genus Tetracanthagyna in these two analyses is not seen with any other genus of the family, and the reason for this discrepancy remains unclear.However, the ITS's sorting of it together with Brachytron fits much better with the morphological and biological criteria [5].
The next four clades each consist of 2-3 specimens of one genus, respectively.These are Rhionaeschna Förster, 1909, represented by three species, Pinheyschna (two species), Anaciaeschna (two species), and Oplonaeschna Selys, 1883 (two specimens of O. armata (Hagen, 1861)).Three of these nodes are well supported, but that of Anaciaeschna is not, and neither did it showed affinity for A. isoceles, sometimes attributed to this genus (see above).
The following robust clade consists of the members of the genus Anax (corresponding to clade 10 in the ITS analysis); however, without A. immaculifrons (see below) and A. ephippiger, both not associated with other Anax and the latter loosely clustering with the genus Triacanthagyna Selys, 1883.The latter strange sorting seems to be an artefact, as it is no longer present if we remove positions in the alignment that are saturated by substitutions or are poorly aligned (Figure S2), and neither is observed in the tree based on the short COI fragment (see below).It was reported before that Anax imperator and the European Anax parthenope s. str.are a rare case of a pair of very different species of Aeshnidae that are not separated based on the COI sequence because of haplotype sharing [10].This partly concerns our COI tree as well (Figure 3), where the sequence 'Anax imperator Germany 16259' does not cluster with the two other sequences of this species (including the specimen from the same place) but appears identical to that of A. parthenope s. str., and so gets to the cluster of the latter.This sequence was obtained by us from a doubtless male specimen of A. imperator.Moreover, the sequence 'Anax imperator Germany 16258' from another male specimen of the same series appeared close to that of A. imperator from Liberia (Figure 3).We have to consider this as a case of introgression of mitochondria from A. parthenope s. str. to A. imperator, similarly to the above-discussed case of the introgression from A. serrata to A. crenata.It is noteworthy that this introgression case was recognised in Germany, whereas Geiger et al. [16] did not register a COI haplotype sharing of A. imperator and A. parthenope in Central Europe.A. imperator and A. parthenope s. str.are clearly separated by the ITS analysis (Figure 2).
Anax parthenope, in the hitherto prevailing broad sense.is not monophyletic in the COI tree, as was also shown before [13], but is split into the West Eurasian and the East Asian branches.The former represents A. parthenope s. str., while the East Asian (including the Far Eastern Russian) specimens represent the taxon A. julius Brauer, 1865.Therefore, we assume the latter as a separate species A. julius, as originally described in detail by Brauer [56] and later again supported by different authors [57][58][59][60].
Strikingly, in the COI tree, Anax immaculifrons forms a lineage which branches from the Aeshnidae stem very early, just after the branching of the Tetracanthagyna clade (Figures 3 and 4).This result appeared to be robust and is reproduced in all our phylogenetic attempts based on the COI gene.This is very strange, not only because it contradicts the well established taxonomy, but also up to the haplotype network (Figure 11), where A. immaculifrons is set apart of other Anax but obviously related to them‚ its root being at the point where Tetracanthagyna branches off.The BLAST search in GenBank [15] for the homology to the COI sequence of A. immaculifrons unequivocally reveals sequences of other species of Anax as most closely related to it.This result is difficult to interpret and the most likely explanation is again an artefact of similarity by chance, which could be facilitated by, e.g., some abnormal substitution rate in the evolutionary lineage leading to A. immaculifrons, or by some structural rearrangement(s).The problem is resolved in the StarBeast analysis (see below).
The last clade is represented by two species and three specimens of Tetracanthagyna waterhousei.In contrast to our ITS phylogenetic reconstruction, it does not cluster with Brachytron, but appears to be the most ancient branch of Aeshnidae.This topology is also robust and is retained in the Gblocks tree reconstructed after the removal of positions by substitutions or due to poor alignment (Figure S1).

Analysis Based on a Short COI Fragment
More sequences are available in GenBank [15] of a shorter (339 bp) fragment of the COI gene, so that more species and genera could be included.
The main topology of the tree reconstructed on its base (Figure 12) did not change as compared to that using the longer COI fragment, with some diverging aspects, as follows.A strange sorting of Coryphaeschna adnexa (Hagen, 1861) in the clade of A. cyanea appeared.But this is no longer present in the Gblock tree (Figure S1), so has to be interpretated as an artefact.The sorting of A. petalura outside the A. cyanea clade seems to be due to shorter sequences of other members in the clade.
In the tree recognised for the shorter COI, some affinity reappeared (but with rather weak support of 0.53) between Caliaeschna, Sarasaeschna, Planaeschna, Cephalaeschna, Remartinia, and Boyeria, so resembling the results of the ITS analysis (clades 13-15, see Figure 2).
The clade including the West Palaearctic genus Brachytron, besides the expected Aeschnophlebia and Epiaeschna, is now updated with Nasiaeschna Selys in Förster, 1900.All these genera share similar biology, being on the wing in spring and preferring strongly vegetated lentic or slowly flowing habitats.For this case, we constructed a haplotype tree with both the short and long fragments of COI, resulting in the same topology.To show also the relationship with the genus Nasiaeschna, the short COI version is given (Figure 13), while the long version is provided in Figure S1.In the short COI fragment tree, Gynacantha africana (Palesot de Belauvois, 1807) and G. vesiculata Karsch, 1891 are found outside of the other Gynacantha, which should be an artifact of insufficient information provided by the short sequence.
In the Anax clade, A. walsinghami, Anax guttatus (Burmeister, 1839), and Anax gibbosulus Rambur, 1842 are now grouped in the first cluster (see above), while the clustering of A. ephippiger with the main Anax cluster has no support (Figure 12).

StarBeast Analysis of COI and ITS Gene Fragments Together
Thus, we have seen that the Bayesian phylogenetic analysis based on the COI fragment using MRBAYES software provided a good resolution on the species level but also some implausible relationships on a higher, genus level.However, it is a matter of fact that any phylogenetic tree based on the investigated gene (the so called 'gene tree') does not necessarily correspond with a phylogenetic tree based on other genes from the same species.Because of insufficient phylogenetic information provided by a particular gene, it is rarely the case that a gene tree is 100% correct.To avoid this problem, StarBeast co-estimates a species tree and several gene trees in one and the same analysis.We have used the StarBeast software to co-estimate species trees based simultaneously on both markers we investigated, COI and ITS.Both entered the analysis for species we sequenced by ourselves, whereas other species entered the analysis with only COI sequences taken from GenBank.We also make species trees where we combined the ITS sequences with the short and the long COI fragments.
The combined analysis by StarBeast of the short COI fragment and ITS region, having more species, as well as that using the long COI fragment and ITS region (Figures 14 and S3) revealed rather a credible topology of the family.The former generally resembles the tree reconstructed with the short COI fragment, as many species entered the joint analysis only with this sequence.In both StarBeast trees, the genus Aeshna (without A. isosceles) is restored as monophyletic, with a good support of 0.83-0.9.These trees also better resolve the clade formed by Rhionaeschna, Anaciaeschna, and Pinheyschna, as already suggested by Ellenrieder [45].Sarasaeschna forms an extra clade away from Planaeschna, Periaeschna, Caliaeschna, and Cephalaeschna, the latter three forming a loose extra clade.Andaeschna is clustered, but without a sound support, in the short COI-ITS tree with Caliaeschna, and forms an isolated branch in the long COI-ITS tree.

Overall Discussion
The main difference between the COI tree and the ITS tree is a greater polytomy of the former, with more basic clades with unresolved phylogenetic relationships, which contain few genera.This could be explained if we suppose that in Aeshnidae, the studied COI fragment probably evolved faster than the ITS region, so that variable positions become saturated with substitutions and the phylogenetic signal is lost at evolutionary distances at which the basic divergence of Aeshnidae took place.
Also, different groups of the genus Aeshna are diverged much deeper than some other well established genera not related to it, e.g., Gynacantha, Heliaeschna, and Austrogynacantha, as well as the genera of Brachytroninae.There are some cases of discordance between the ITS and COI trees.The most striking is the position of the genus Tetracanthagyna, which in the ITS tree is among other Brachytroninae (Figure 2), as also suggested by the morphology-based phylogenetic analysis [5], but appears to be the most ancient divergence of the Aeshnidae stem on the COI tree (Figure 3).This discordance also persisted in the joint StarBeast analysis (Figures 13 and S3).The second discordance concerns Austrogynacantha, which was placed in the ITS based tree as expected, with Gynacantha (Figure 2), but had no relation to Gynacantha in the COI-based trees.The next evident problem in the COI tree was the position of A. immaculifrons, which clusters outside Anax in the COI tree (Figure 3).The joint StarBeast analysis resolved this problem and reintegrated this species into Anax (Figures 5 and S3).
At the low taxonomic level of species, we may point to the discordant phylogenetic pattern in the A. juncea/A.subarctica group, wherein the ITS sequences separated these two species (Figure 2), while the COI sequences also formed two clades, but one of them included both A. subarctica and A. juncea from eastern Eurasia and North America and the other included A. juncea from western Eurasia (Figure 3). A. parthenope does not diverge from A. julius in the ITS tree (Figure 2) or from A. imperator in the COI tree (Figure 3), but not vice versa.
Extensive molecular phylogeny of insects at the levels of species and genera started using the mitochondrial gene COI as the most popular marker, which was suggested for insect barcoding.It had such advantages as fast evolution, sometimes allowing researchers to trace divergence even at intra-species level (e.g., A. cyanea and A. juncea), existing in numerous copies of the cell and strictly (in animals, with few exceptions) maternal inheritance excluding recombination.With time, evidence accumulated that mitochondria evolution may be oddly discordant to that of the nuclear genome and hence phenotype [61].This concerns Odonata as well, with most examples coming from Coenagrionidae [49,50].In our analysis we found a possible case of introgression between A. serrata and A. crenata, evidenced by two specimens, from Finland and the Ural Mountains, which were morphologically A. crenata and had the ITS sequences of A. crenata, but clustered in the COI analysis within A. serrata (Figure 3).Another putative case of introgression of COI was supposed to take place from A. parthenope to A. imperator.

Taxonomic Implications
Our phylogenetic trees based on ITS sequences suggest that the genus Aeshna in the current sense is not monophyletic.According to the ITS tree (Figure 2), it can be made monophyletic if we synonymize with it the genera Pinheyschna and Polycanthagyna (the next node of which unites the current Aeshna spp., and Pinheyschna only has a weak support of 0.68).Even this broader solution would still place A. isoceles outside Aeshna.At the same time, in both COI trees (Figures 3 and S2), Aeshna is monophyletic (although with a weak support), but again without A. isoceles.In none of our analysis did A. isoceles show a closer relationship to the genus Anaciaescha, as has sometimes been suggested before [43,55].Therefore, we had to erect a new genus solely for A. isoceles.This is not surprising, as this has been discussed for about 100 years, when Friedrich Ris asked Erich Schmidt "What is Aeshna isosceles?"[43].However, the placement in Anaciaeschna, as suggested by him [43], cannot be followed, as none of our gene fragments investigated by different algorithms correspond to this assumption.Our results (Figures 2 and 3) also did not support a relation of A. isosceles to the genus Andaeschna De Marmels, 1994, as discussed by von Ellenrieder [45].
We, therefore, suggest for A. isoceles the new genus: Isoaeschna gen.nov.Type species: Libellula isoceles Müller, 1767.Ethymology: 'ίσoς'-a Greek prefix meaning 'equal', Aeshna-the name of the genus to which the type species was attributed for a long time.
Differential diagnosis (based on [43,45]): This monotypic genus has some unique combinations of morphological and colourational characters: a transverse ridge on the sternum of S2 (shared with Anaciaeschna); narrow, parallel-sided auriculae (shared with Anaciaeschna); rounded hindwing anal angle (shared with Andaeschna and Anax); anterior and posterior veins of anal triangle fused at a point, without prolongation of the fused vein (shared with Andaeschna and Anax); membranule length comprising 75-100% of the wing anal margin (shared with Anaciaeschna and Andaeschna) [45]; green eyes without any trace of blue, presence of amber hindwing basal spots [43], the absence of the T-marking on the frons (present in most genera of Aeshnidae but also absent in Andaeschna and most Anax) [45]; a yellow dorsal triangle on S3.According to von Ellenrieder [45], the new genus appears most close to Andaeschna, differing from it in the presence of the transversal ridge on the S2 sternum (a conical tubercle bearing denticles in Andaeschna) and the parallelsided auriculae (triangular or quadrangular with denticles in Andaeschna).However, our molecular analysis did not prove this affinity.
Polycanthagyna and Indaeschna cluster together very closely in the COI trees and in the combined gene analysis by StarBeast.The ITS sequence of Indaeschna is still missing; but we nevertheless suggest to synonymize these two genera: Indaeschna Fraser, 1926 = Polycanthagyna Fraser, 1933, syn.nov.All our results are unequivocally in favour of subsuming the genera Aeschnophlebia Selys, 1883, Nasiaeschna Selys in Förster, 1900, and Epiaeschna Hagen in Selys, 1883 under the genus Brachytron Evans, 1845, which has so far been monotypic.This is also not surprising and was already discussed by others [5,51].All these dragonflies show great similarity in morphology and autecology.Thus, the following synonymies are put forward: Brachytron Evans, 1845, valid name = Aeschnophlebia Selys, 1883 syn.nov.= Epiaeschna Hagen in Selys, 1883 syn.nov.= Nasiaeschna Selys in Förster, 1900 syn.nov.
More complicated and unresolved remains the position of the genus Tetracanthagyna.While in the ITS tree, it clusters with Brachytron, as discussed earlier [5]; in all COI analyses, it was sorted outside as a primeval clade.
Caliaeschna clusters with Periaeschna and Cephalaeschna in the ITS tree (Figure 2), but in the long COI tree the clustering of Caliaeschna with other genera is too loose and the closest genus is Sarasaeschna (Figure 3), for which we have no ITS sequence.Thus, our results are too equivocal for a definitive taxonomic merger of Caliaeschna and Cephalaeschna, as suggested earlier [5,48].
The genus Anax was found to be monophyletic in the ITS tree (Figure 2), whereas in the COI trees (Figures 3 and 12), two species fall outside: A. ephippiger is not clustered or loosely clusters with the rest of Anax, while Anax immaculifroms is found near the base of the tree.However, the joint StarBeast analysis restores its position among Anax.
Boyeria unequivocally forms an extra clade in all our analyses.In some genera, a deeper divergence can be recognized, so, for example, the African and Asian members of the genus Gynacantha form two subclades in the COI trees, respectively.
The genera Anaciaeschna, Rhionaeschna, and Pinheyschna are in the same clade in the COI and StarBeast analysis and may be regarded as related, despite their geographical separation, as discussed earlier [58].Unfortunately, we had no sequences of the genera Zosteraeschna Peters at Theischinger, 2011 to check to see if they would belong to the same clade, as expected.
Some taxonomic inferences at the species level can be made.Anax julius is well separated from A. parthenope in the COI tree (Figure 3) but not in the ITS tree (Figure 2); the former is in favour of there being different species, as proposed earlier [56][57][58][59][60]62].The lack of divergence of their ITS region could be ascribed to the above-mentioned putative slower evolution of the ITS region in Aeshnidae.
All our analyses unequivocally suggest that the North American A. septentrionalis and the Eurasian A. caerulea are extremely close to each other.As has been repeatedly stated [40,47,48,[62][63][64], they have no substantial morphological differences, while the reported ones were scarcely distinctive.Therefore, following [47,48], we treat the American populations as the subspecies Aeshna caerulea septentrionalis.
The American species A. palmata and A. constricta share the most common COI allele (Figure 14) and look like the same species.We, however, abstain from their synonymization, as both are known to broadly co-occur and to differ by a number of characters.Maybe our result reflects some mitochondrial introgression between these species.
We synonymize Gynacantha hyalina Selys, 1882 syn.nov.with Gynacantha subinterrupta Rambur, 1842.G. dravida Lieftinck, 1960 looks like the same species as Gynacantha subinterrupta Rambur, 1842 in the COI and haplotype tree; however, we did not have the ITS sequence to decide this definitively.
The two recently proposed [25], closely related but separate species, A mixta and A. soneharai, are well separated in the ITS tree (although by three substitutions only) (Figure 2) and the COI haplotype networks (Figure 8), while in the COI tree the latter looks like an in-group inside the former.
The ITS analysis of A. juncea revealed a separation of A. juncea from the Caucasian/Transcaucasia region versus elsewhere.This may deserve taxonomical fixation at the subspecies level.Two available names were proposed for A. juncea from the Caucasus: A. juncea crenatoides Bartenev, 1925 and A. juncea atshischgho Bartenev, 1929 [65-67].They were claimed to share such unfortunately quantitative characters as broad thoracic stripes and shallowly incised vulvar lamina, and to differ in the absence (in the former) versus presence (in the latter) of the so-called 'lateral genital plates' in the female ovipositor [67].It is noteworthy that we managed to involve into our analysis specimens from a population where those 'lateral genital plates' were present, from North Caucasus, and from populations where they are absent, from Dagestan [68] and Georgia (from where A. juncea crenatoides was described) [52].Both their ITS and COI sequences appeared identical.Although subspecies are entities of geographical variation usually differing in some single character and so do not need to be diverged all over their genomes, our result is in favour of treating these subspecies as synonyms; A. juncea crenatoides = A. juncea atshischgho syn.nov.
The divergence of the COI gene of the same A. juncea by longitude, from the West Palaearctic east to West Siberia and America plus the East Palaearctic west to Pakistan, with the species A. subarctica clustering to the latter, as can be seen in Figure 3, is striking.In the ITS tree, A. juncea and A. subarctica perform as well diverged monophyletic species (Figure 2).We may suppose the following scenario which could have taken place during the repeated coolings and warmings of the Pleistocene/Holocene. First, both species diverged from their common, most probably Eurasian ancestor in different continents to become A. juncea in Eurasia and A. subarctica in North America.Then, after some of the repeated restorations of Beringia, both expanded to the other continent.The expansion of A. juncea to America was accompanied by mitochondrial introgression from A. subarctica to A. juncea.Then, those 'contaminated' populations of A. juncea expanded back to Eurasia, occupying its eastern regions.This is more or less concordant with the results of an attempt at phylogeographical analysis of the same COI gene by the same two species [11], but the data lost most geographical information due to operating in such a huge 'region' as 'Russia', which occupies more than half of the Holarctic, so they hardly provide an informative geographical resolution.
Our suggestions for taxonomical changes in Aeshnidae, as discussed above, are summarized in Table 2.

Figure 2 .
Figure 2. Bayesian phylogenetic tree reconstructed from the ITS region of representatives of Aeshnidae using MRBAYES 3.2.7a,shown as a phylogram on the left side and the mirrored cladogram on the right side.Bayesian posterior probability values are depicted at the nodes.Included are our own sequences (PCR number next to the name) and those retrieved from GenBank (accession number next to the name).Clade 3 (cyanea-clade) includes Aeshna cyanea (Müller, 1764), Aeshna vercanica Schneider et al., 2015, and the East Asian Aeshna petalura Martin, 1908.Males of all these three species have a particular morphology of the upper appendages, which are broad and have a downward terminal hook, which resembles a raptor's beak.The possible relationship of A. petalura with A. cyanea was already suggested by Erich Schmidt [43].Two clusters of A. cyanea can be recognized, one including specimens from the Caucasus and Transcaucasia and the other with specimens from Europe and North Africa (Figure2, see the cladogram on the right side).This dichotomy has been recognized before[14,44], and is paralleled by the above-discussed divergence of A. juncea from the Caucasus and Transcaucasia versus from elsewhere.The American members of the 'umbrosa group' ('paddle-tails') were not available for the ITS analysis.Clade 4 includes Aeshna mixta Latreille, 1805 and Aeshna soneharai (Asahina, 1988) as very closely related sister taxa.Actually, their ITS region differs in just two substitutions in the ITS region.It is noteworthy that Onishko et al.[23], who raised the taxon A. soneharai to species level based on external characters, behaviour, sympatric occurrence, and differences in the mitochondrial COII sequences, also sequenced the ITS2 spacer and found it to be identical to A. mixta.We now have sequenced the broader ITS region and found that Aeshna soneharai differs in three nucleotide substitutions from A. mixta.Clade 5 includes Aeshna affinis Vander Linden, 1820 only.Clade 6 is represented by the two species of the genus Anaciaeschna Selys, 1878, Andaeschna occidentalis Bota-Sierra, 2019, and Pinheyschna yemenensis(Waterston, 1985).This clade represents members of three continents.Such a relation was already discussed by von Ellenrieder[45].Clade 7 contains both species of the genus Polycanthagyna Fraser, 1933 available.Clade 8 is represented solely by the species Aeshna caerulea (Ström, 1783).

Figure 3 .
Figure 3. Bayesian tree reconstructed from the long (632 bp) fragment of the COI gene of representatives of Aeshnidae using MRBAYESs 3.2.7a,shown as a phylogram on the left side and the mirrored cladogram on the right side.Bayesian posterior probability values are depicted at the nodes.Included are our own sequences (PCR number next to the name) and those retrieved from GenBank (accession number next to the name).

Figure 4 .
Figure 4. Haplotype network of the long COI fragment for Aeshna cyanea, A. vecranica, C. umbrosa, A. constricta, A. palmata, and A. petal.Like in the ITS analysis, specimens of A. cyanea from North Africa and Europe cluster together.However, specimens of A. cyanea from the Caucasus and Transcaucasia region on one hand and from the rest of the range on the other had no longer form sister clades as in the ITS tree (Figure2).Instead, the Caucasian plus Transcaucasian specimens radiate from the base of the A. cyanea cluster, while the European plus African cluster is now internal, as a sprouting among them.This can be interpreted as the species A. cyanea having originated in the Caucasian area and then one of its lineages having spread to the west and occupied vast European and North African territories.A. vercanica and A. cyanea are sister branches in the COI tree (Figure3), but the node of A. cyanea has a weak support of 0.55.

Figure 5 .
Figure 5. Haplotype network of the long COI fragment for Aeshna caerulea, A. septentrionalis, and A. sitchensis.

Figure 6 .
Figure 6.Haplotype network of the long COI fragment for Aeshna juncea and A. subarctica.

Figure 9 .
Figure 9. Haplotype network of the long COI fragment for representatives of the genus Gynacantha.
Two clusters can be recognized in the main Anax clade: one consists of Anax tristis Hagen, 1867, A. junius (Drury, 1773), A. julius, A. nigrofasciatus Oguma, 1915, A. imperator and A. parthenope; the second consists of A. congoliath Fraser, 1953, A. gladiator Dijkstra et Kipping, 2015, A. speratus Hagen, 1867 and A. rutherfordi McLachlan, 1883.The same relationships between the Anax species are seen in the haplotype network (Figure 11).In the combined gene analysis by StarBeast (see below), Anax walsinghami McLachlan, 1883 is added to the first group, and A. immaculifrons is sorted between all Anax species and A. ephippiger.

Figure 11 .
Figure 11.Haplotype network of representatives of the genera Anax, Triacanthagyna, and Tetracanthagyna for the COI long fragment.

Figure 12 .
Figure 12.Bayesian tree reconstructed from the short (341 bp) COI gene fragment of representatives of Aeshnidae using MRBAYESs 3.2.7a.Bayesian posterior probability values are depicted at the nodes.Included are our own sequences (PCR number next to the name) and those retrieved from GenBank (accession number next to the name).

Figure 14 .
Figure 14.Multi-locus sequence species tree reconstructed with StarBeast3 v 1.1.7based on the short COI gene fragment and the ITS region of representatives of Aeshnidae.Bayesian posterior probabilities values are depicted at the nodes and as colour in the branches.

Table 1 .
Information on Aeshnidae specimens used; the boldfaced GenBank reference numbers refer to sequences obtained in the course of this study.

Table 2 .
Valid names according to the taxonomic treatments of the present work.