Since 1999, the entomologists and orthopterologists Baudewijn Odé, Rob Felix, Luc Willemse, and Roy Kleukers have been organizing Orthoptera collection expeditions on the field, in southern Europe, focusing primarily on Spain, Portugal, Italy, Greece, and Romania. Voucher specimens were killed using ethyl acetate, dried, mounted, and identified utilizing literature and identification keys specifically designed for European orthopterans (e.g., Willemse et al. 2018; Iorio et al. 2019). Some individuals were included in the dataset even if they could not be identified to species level. Based on the barcode data results, more precise information will be possible in future projects. The right middle legs were taken from freshly killed specimens and preserved in 97% ethanol (Rentz 2010; Willemse et al. 2018). A total of 1,698 voucher specimens are stored at the collections at the Naturalis Biodiversity Center, the national research institute for biodiversity located in Leiden, Netherlands.

The DNA extraction, amplification and sequencing were performed in two parts. 630 samples were processed at the Leibniz Institute for the Analysis of Biodiversity Change (LIB) in Hamburg (Germany), and 1,068 samples at the Naturalis Biodiversity Center in Leiden (Netherlands).

At LIB, DNA was extracted according to the methods of Chelex® 100 resin-based protocol (de Lamballerie et al. 1992). For PCR amplification of the DNA barcode region, the COBL (forward) and COBU (reverse) primers were used, tailored for Orthoptera to improve amplification and sequencing results (Huang et al. 2013). The reaction master mix composition and cycling protocol are outlined as supplementary material (File S1 [protocol a]). The PCR products were checked for successful amplification on an agarose gel and recorded using a gel imaging system. Samples were sequenced by Macrogen Europe (Amsterdam).

The Sanger sequencing results (ab1 files) underwent a quality check. Sequences chromatograms were analysed using Geneious Prime 2023.1.1, trimming poor-quality edges and manually correcting bases with low-quality peak maps. Sequences with unreliable peaks were subjected to PCR and sequencing repetition, adjusting the annealing temperature to 49°C, which improved outcomes in several instances. Sanger sequencing data were filtered, selecting sequences longer than 500 bp (base pairs), except for some species where only few longer sequences were available.

At Naturalis, DNA extraction was performed adhering to the ARISE (Authoritative and Rapid Identification System for Essential biodiversity information) protocol (van Ommen Kloeke 2022; protocols.io; bomb.bio). Following DNA extraction and cleanup using KingFisher, PCR was conducted on all samples. This step employed the same primers (COBL and COBU).

PCR products were verified on a 2% agarose E-Gel, with nearly all samples showing positive results, prompting continuation with all samples (detailed PCR protocol in File S1 [protocol b]). Subsequently, the Nanopore dual barcoding protocol with kit 14 was applied for the sequencing. In this protocol, samples with no bands in the elektrophoresis were also used for library preparation.

All sequences were aligned using MAFFT V 7.505 (Katoh and Standley 2013) on the CIPRES Science Gateway (Miller et al. 2010). Sequences with excessively wide gaps and those containing internal stop codons were removed. Next, all sequences underwent contamination or misidentification checks using BLASTn (Camacho et al. 2009).

The newly generated sequence data and metadata were then uploaded to the Barcode of Life Data Systems (BOLD PROJECT “MEDOR Barcoding of Mediterranean Orthoptera”). Thanks to our preliminary quality assessments, all barcode sequences met the automatic quality criteria of BOLD upon upload.

The initial dataset was augmented with two additional datasets of European Orthoptera sequences available on the BOLD system, named GBORT-GBOL (Hawlitschek et al. 2016; 745 records) and IBIOR (Pina et al. 2024; 420 records), resulting in the merged dataset DS-MEDOR1. The former consists of records primarily from Austria, Germany, and Switzerland, while the latter comprises records from continental Portugal. This integration aimed to achieve a more accurate subdivision of the samples by species within the phylogenetic tree. This was facilitated by the fact that the datasets from other European and Mediterranean basin areas share several species and genera with those in the current study.

DNA barcoding shows variable efficiency between Caelifera and Ensifera due to factors such as numts frequency (Kaya and Çıplak 2018), hybridization rate (Hawlitschek et al. 2016), ILS (Nabholz 2023), and Wolbachia infections (Ilinsky et al. 2022). Thus, phylogenetic and species delimitation analyses were conducted separately for each suborder (unless stated otherwise), which are well-supported as monophyletic lineages (Song et al. 2015; Cigliano et al. 2024), dividing the merged dataset into two subsets. This division significantly reduced the initial dataset size, facilitating the management of the extensive number of sequences.

Two preliminary Maximum Likelihood (ML) trees (Felsenstein 1981) were reconstructed of the aligned dataset using IQ-Tree v2.2.2.7 (Nguyen et al. 2015). The GTR (General Time Reversible) + I + G model (Tavaré 1986), identified by ModelFinder (Kalyaanamoorthy et al. 2017) as the best-fit model for nucleotide substitution, was used in both analyses. The merged datasets revealed sequences occupying potentially incorrect or ambiguous positions within the phylogenetic trees. To resolve possible issues of contamination or misidentification, voucher specimens were re-verified with the assistance of experts at the Naturalis Biodiversity Center in Leiden. In accordance with the approaches recommended by Song et al. (2008) and Leite (2012), chromatograms and sequence editing were meticulously re-examinated. Furthermore, bioacoustic, biogeographic, and morphological data were consulted to ensure accurate species verification. Sequences presumed to represent numts or contaminations were excluded from the final datasets and trees.

Specimens not identified at the species level were retained for further analysis, provided their most specifically identified clade was accurately placed within the ML tree, aligning with currently accepted Orthoptera taxonomy (Cigliano et al. 2024).

After quality measures, each definitive subset was re-aligned and trees reconstructed using the same parameters as for the preliminary analyses described above.

All phylogenetic analyses were executed with XSEDE (eXtreme Science and Engineering Discovery Environment) through the CIPRES Science Gateway (www.phylo.org; Miller et al. 2010).

The effectiveness of DNA barcoding in species identification was assessed by calculating the ratio between the number of species not exhibiting BIN sharing and the total number of species assigned a BIN code, following the methodology described by Hawlitschek et al. (2016). Additionally, a weighted average was computed to account for the varying number of samples per genus, ensuring that genera with a larger number of samples had a proportional influence on the final estimate (Sokal and Rohlf 1995). The weighted average of species per genus was calculated by multiplying the number of species within each genus by the number of samples representing that genus, summing these products, and dividing by the total number of samples across all genera. Both statistics were calculated for the entire dataset and for caeliferans and ensiferans separately.

Species delimitation analyses were conducted on the definitive datasets using various species delimitation approaches to detect and compare MOTUs, to explore species diversity and deepen the understanding of species boundary delimitation issues. For this purpose, both similarity (BIN, ABGD, ASAP) and clustering-based (GMYC, PTP) approaches designed for a single-locus strategy were utilized. In all methods we employed, the term “partition” refers to the grouping of sequences into distinct clusters or units, typically representing putative species. In methods such as ABGD and ASAP, partitions emerge from the identification of genetic gaps or distance thresholds, while in GMYC and PTP, partitions are determined by tree-based criteria, such as branch length distributions or coalescent processes (Pons et al. 2006; Puillandre et al. 2011, 2020; Ratnasingham and Hebert 2013; Zhang et al. 2013).

BIN method: the Cluster Sequences tool, implemented by BOLD, generates MOTUs through the REfined Single Linkage (RESL) algorithm, grounded on uncorrected pairwise distances (p-distance). The Barcode Index Number (BIN) informatics system provides a unique alphanumeric code for MOTUs (BOLD: 3 letters, 4 numbers (Ratnasingham and Hebert 2013)). BIN-assignment data were retrieved from the dedicated web interface for all chosen specimens. Cases of BIN discordance arise when traditionally recognized species encompass more than one BIN, whereas BIN sharing pertains to single BINs containing members from multiple recognized species (Hebert et al. 2004).

ABGD method: the Automatic Barcode Gap Discovery was adopted to divide samples based on genetic distance, detecting the so-called “barcode gap” (Puillandre et al. 2011). The analysis used a relative gap width of 1 and 20 bins, and results were interpreted using a prior intraspecific divergence limit of P = 0.01 (Puillandre et al. 2011; Gonçalves et al. 2021) with four combinations of different substitution models and partitions employed: Jukes-Cantor (Jukes and Cantor 1989) or Kimura two-parameter (Kimura 1980) models and initial or recursive partition.

ASAP method: the Assemble Species by Automatic Partitioning creates new partitions by amalgamating sequences at equal pairwise distances into progressively larger groups until the final partition encompasses all records (Puillandre et al. 2020). This analysis used Jukes-Cantor distances and default parameters, with the best initial and recursive partitions selected based on the lowest “asap-score”.

For both ABGD and ASAP, the initial partition refers to the first grouping of sequences based on broader genetic thresholds, while the recursive partition further refines these clusters by progressively splitting them to detect finer genetic distinctions, potentially identifying cryptic species. Both ABGD and ASAP analyses were run on the respective online platforms (available at https://bioinfo.mnhn.fr/abi/public/abgd and https://bioinfo.mnhn.fr/abi/public/asap/asapweb.html, respectively).

GMYC method: the Generalized Mixed Yule Coalescent discerns the most likely putative species clusters basing on an ultrametric tree (Pons et al. 2006). Concatenated sequence alignments were analysed using BEAST v.2.6.7 (Drummond and Rambaut 2007; Suchard and Rambaut 2009; Bouckaert et al. 2014), a Bayesian inference tool generating ultrametric trees. Parameters included a range of clock and speciation model priors (strict or optimized relaxed clock model and Yule or Coalescent Constant Population process). According to ModelFinder (implemented by IQ-Tree), GTR+I+G was chosen as the most suitable model of nucleotide substitution. Due to predominantly non-converging runs and lower-than-expected Effective Sample Size (ESS) values, likely from over-parametrization, the partition scheme was adjusted to employ the simpler HKY model (Hasegawa et al. 1985). Each run, totaling eight, was set for 100 million iterations, adjusting the burn-in rate for optimal convergence. Tree and parameter sampling occurred every 1,000 steps, with TRACER v.1.7.2 employed to ensure ESS values exceeded 200, as advised by Drummond et al. (2007). The posterior distribution was summarized into the maximum clade credibility tree using TreeAnnotator v.2.6.6 (Drummond and Rambaut 2007). Newick files, post-burnin, were then subsampled to 10,000 trees each for processing. The resulting ultrametric trees supported only single-threshold GMYC analyses to prevent species number overestimation, a common issue with multiple-threshold approaches (Fujisawa and Barraclough 2013). The GMYC species delimitation analysis was then executed within R Software (Team, R Core 2014), using the “ape” package for phylogenetic analysis (Paradis and Schliep 2019), “paran” for parallel analysis (Dinno 2012), “rncl” for reading and manipulating phylogenetic data (Michonneau et al. 2016), and “splits” for species delimitation and network analysis (FitzJohn et al. 2009).

PTP method: this method delineates hypothetical species clusters by analysing branch length distributions in a rooted non-ultrametric gene tree. The PTP model, enhanced by the bPTP version with Bayesian support for delimited species on the input tree (Zhang et al. 2013), underwent analysis via an online server for both PTP and bPTP versions (http://species.h-its.org), targeting the maximum likelihood tree as input. Analysis parameters included 500,000 MCMC generations, a thinning of 500, and a burn-in of 0.1, with convergence checked as Zhang et al. (2013) recommended.

The Python-based script SPdel v.2.0 (Ramirez et al. 2023) was utilized to visualize a summary representation of all delimitation methods, including the ML tree, specimen labels, and putative species clusters. This was computed according to each of the five species delimitation methods implemented and their different combinations of parameters, along with a manually designed consensus, which compares features summarizing all results based on a majority criterion (MOTUs found in more than 50% of methods). Records lacking a BIN code were excluded from the species delimitation summary, as SPdel can only provide graphical output for records assessed by all involved methods. All species were still represented except for Glyphanus obtusus Fieber, 1853, Tettigonia caudata (Charpentier, 1845), and Isophya lemnotica Werner, 1932. For these species, the delimitation analysis was also conducted, but results were only described and not displayed in the summary representation.

A comprehensive description of the command, input files, and options used for the SPdel.py analysis can be found in File S8.

The final MEDOR (Barcoding of Mediterranean Orthoptera) alignment comprises 1,441 newly generated, quality-checked barcodes, each 726 bp in length (accession numbers for newly generated sequences on BOLD (Ratnasingham and Hebert 2013) are listed in File S5). These barcodes correspond to 270 identified species, including 10 subspecies. Additionally, 209 barcodes lack species-level identification, with most assigned to the genus level (49 genera). Among these, three records are identified as Pamphagidae sp., two as Ephippigerini sp., and one as Tettigoniinae sp. Notably, the dataset includes barcodes for 26 species (and two subspecies) of ensiferans and four species (and one subspecies) of caeliferans that had no publicly available barcodes before this study.

Incorporating data from Hawlitschek et al. (2016) and Pina et al. (2024), the complete dataset totals 2,606 specimens across 351 identified species, with a total of 14 subspecies and an equivalent number of taxa not identified at the species level (209). A significant portion of specimens were collected in Spain and the Canary Islands (851), Germany (544), Portugal (511), Italy (173) and Greece (145), followed by 15 other countries in the Mediterranean Area, plus two specimens from Sri Lanka and two from the United Arab Emirates (Fig. 1, with the exclusion of non-European specimens). Relative to the known Orthoptera species richness described in each country (Cigliano et al. 2024), the coverage of Orthoptera biodiversity is 100.00% in Germany, 93.44% in Portugal, 63.47% in Spain, 43.89% in Italy, and 31.12% in Greece (further details in File S5).

Figure 1. 

Sampling localities of all European samples involved in this study. Color legend: blue [MEDOR project], green [GBORT-GBOL project (Hawlitschek et al. 2016)], orange [IBIOR project (Pina et al. 2024)]. The four samples from Sri Lanka and Saudi Arabia are not included in this map.

The analysis of the merged dataset revealed significant variation in barcoding effectiveness between ensiferans and caeliferans. Overall, 76.72% of all studied Orthoptera species were accurately identified through DNA barcoding. However, the success rate differs markedly between the suborders: 63.64% for Caelifera, primarily due to BIN sharing among certain species, compared to 88.53% for Ensifera. The datasets for ensiferans and caeliferans also highlight differences in the species richness of the respective genera, with a weighted mean of 3.78 and 8.29 species per genus, respectively.

The ML trees of Caelifera and Ensifera are shown in Figs 2 and 3, and more in detail in File S2 [graphs a, b] respectively. Species delimitation analyses were then performed separately on two distinct datasets: Caelifera, with 1,478 records, and Ensifera, with 1,128 records.

Figure 2. 

Maximum Likelihood tree of Caelifera from Central and Southern Europe. The tree was created with IQ-Tree v.2.2.2.7 on XSEDE (eXtreme Science and Engineering Discovery Environment) through the CIPRES Science Gateway (www.phylo.org). The tree displays 1,443 barcodes of identified (164) and unidentified species. Higher taxonomic levels (all families and the subfamilies of Acrididae) are coloured in the tree. The species depicted inside the tree is Arcyptera (Arcyptera) tornosi Bolívar, 1884.

Figure 3. 

Maximum Likelihood tree of Ensifera from Central and Southern Europe. It was created with IQ-Tree v.2.2.2.7 on XSEDE (eXtreme Science and Engineering Discovery Environment) through the CIPRES Science Gateway (www.phylo.org). The tree displays 1,111 barcodes of identified (185) and unidentified species. Higher taxonomic levels (all families and the subfamilies of Tettigoniidae) are coloured in the tree. The species depicted inside the tree is Pseudomogoplistes vicentae Gorochov, 1996.

The various species delimitation approaches employed yielded hypothetical OTU counts ranging from 355 to 505, compared to 351 previously identified species (Table 1). BOLD identified 445 BINs (162 for caeliferans and 283 for ensiferans). 278 species demonstrated BIN concordance, where all barcodes from specimens of a given species grouped into a single BIN, exclusively containing barcodes from that species alone. 95 species displayed singleton BINs, characterized by containing only a single record without additional sequences. BIN discordance was observed in several species, indicating a mismatch between the BIN system and traditional taxonomic classifications. 21 and 10 BINs, in caeliferans and ensiferans respectively, were identified as shared across multiple species, impacting 60 and 21 species in total. Additionally, 35 and 17 records were not assigned any BIN due to the absence of sequences longer than 500 bp, the minimum length required for BIN assignment on BOLD. These records were omitted from the BIN discordance analyses. However, for each of these, there was at least one co-specific record that did possess a BIN code, with the sole exceptions being Glyphanus obtusus, Isophya lemnotica, and Tettigonia caudata. Detailed associations between BIN codes and the associated species are provided in File S4 [tables a, b] and in File S6.

The ABGD method, especially its initial partition using the Kimura two-parameter model, emerged as the most conservative, contrasting with the PTP method, which identified the highest number of hypothetical species. Parameter variations within ABGD resulted in partitions of varying sizes, with the Jukes-Cantor model generally less conservative than the Kimura two-parameter, leading to a higher MOTUs count across both suborders, apart from the initial partition for caeliferans. Recursive partitions consistently indicated more species than initial partitions. Contrary to ABGD findings, the ASAP method with Kimura two-parameter resulted in larger partitions compared to those from Jukes-Cantor.

The GMYC method, applied with four different parameter combinations, produced a consistent number of species clusters for both ensiferans and caeliferans. However, the application of the Coalescent Constant Population process combined with a relaxed clock model resulted in an anomalously high species count for this dataset, leading to its results being set aside due to poor fit. The outcome of this specific combination, indicating an unexpectedly large number of species, suggests it may not be suitable for this particular dataset analysis. The PTP output in caeliferans surpassed the average species count across BIN, ABGD, ASAP, and GMYC methods, while in ensiferans, the 271 clusters of PTP aligned closely with the averages from other methods. The bPTP analysis, despite extensive generation counts, failed to converge, leading to the exclusion of its results in favour of maximum likelihood outcomes, as recommended by Zhang et al. (2013).

The graphic outputs from SPdel are depicted in File S3 [graphs a, b], showcasing delimitation bars in the order of 11 method and parameter combinations, culminating in a final consensus.

In the maximum likelihood tree, while barcodes are correctly assigned in the absence of BIN sharing, higher taxonomic groups are often not retrieved as monophyletic. This suborder’s species delimitation methods generally concur on the number of putative species clusters, except for PTP, which delineates a higher species count. The hypothetical species range from 118 to 251. For further details see Table 2.

In the maximum likelihood tree, while barcodes are correctly assigned in the absence of BIN sharing, higher taxonomic groups are often not retrieved as monophyletic. In Ensifera, the different delimitation methods show an overall agreement on the number of putative species clusters, with the only exception of PTP, showing a higher number of delimited species. The hypothetical species range from 118 to 251. For further details see Table 3.

The percentages of DNA barcoding effectiveness in species identification are somewhat lower than those reported by Hawlitschek et al. (2016) (78.2% for Caelifera + Ensifera, 59.1% for Caelifera, and 100% for Ensifera), but the current findings are based on a substantially larger dataset. Furthermore, in Hawlitschek et al. (2016), each ensiferan genus was often represented by only one species, making it less likely to include any closely related species that might share barcodes. The datasets for ensiferans and caeliferans also highlight differences in the species richness of the respective genera. This significant difference suggests that the higher species richness within each caeliferan genus could contribute to the lower barcoding effectiveness observed, as genera with many species are more likely to exhibit BIN sharing and overlapping barcode clusters.

As expected from a barcode tree (Vences et al. 2005; Hawlitschek et al. 2016), the inferred ML tree showed the taxonomic relations above the generic level as poorly adherent to the recognized taxonomy.

The detailed discussion of the delimited putative species clusters is organized according to traditional taxonomic families and can be found in File S7 [part 1] for caeliferans and File S7 [part 2] for ensiferans.

The employment of multiple species delimitation methods, each with distinct parameter combinations, on a comprehensive dataset has facilitated the identification of various patterns of diversity, some of which are clearly delineated, while others remain unresolved. In total, 12 cases for caeliferans and 42 cases for ensiferans were identified as potentially harboring more than one cryptic lineage. Additionally, there were 10 cases for caeliferans and 26 for ensiferans of unidentified independent MOTUs. The detected instances of barcode sharing (12 cases for caeliferans and nine cases for ensiferans) can either be attributed to intraspecific variability, suggesting a need for synonymization among these taxa, or may highlight the ineffectiveness of using a single genetic marker for accurate species delimitation, like the case of Chorthippus spp.

Consistent with prior studies, delimitation patterns within the Gomphocerinae (Hafayed et al. 2023) and Oedipodinae (Kock et al. 2024) subfamilies (Acrididae) are predominantly influenced by recent radiation events. Species emerging from such events often evade detection through DNA barcoding, resulting in shared BINs and indistinguishable barcodes across all delimitation algorithms. Notably, analogous patterns are also evident in some Tettigoniinae, particularly within the genera Platycleis Fieber, 1853 and Tessellana Zeuner, 1941 (Platycleidini), which are believed to have undergone recent divergence (Mugleston et al. 2018).

Among well-represented taxa of Caelifera, instances of potential cryptic diversity are notably prevalent within the Pamphagidae and Tetrigidae families, as well as the Calliptaminae and Melanoplinae subfamilies of Acrididae. More distinct and frequent occurrences of potential cryptic diversity are observed within Ensifera, especially among Gryllidae, and also Tettigoniinae, Bradyporinae, and Phaneropterinae (Tettigoniidae). Ensiferans encompass a greater number of taxonomically recognized species (185 vs. 164) and the highest number of putative species clusters identified by each delimitation method, despite a smaller initial specimen count (1128 vs. 1478).

Among all the delimitation methods compared, ABGD, particularly when using the combination of Kimura two-parameters and the initial partition, produced results most consistent with traditional taxonomy. On the other hand, the PTP method was the least consistent, oversplitting the lineages into a large number of clusters.

Geographical regions partially reflect genetic diversification. Various potential cryptic taxa have indeed been identified in isolated localities, such as islands and glacial refugia. Key examples are represented by Anterastes supersp. serbicus Brunner von Wattenwyl, 1882 in Anatolia and the Balkans (Çiplak et al. 2010a), but also the genus Oedipoda Latreille, 1829 in the Mediterranean (Hochkirch et al. 2023), both showing correlation between the radiation of the group and the topography of the species’ ranges. The diversification and frequency of these taxa are intimately linked to climatic and geological events, such as Pleistocene glacial cycles, which fragmented populations and created ecological niches for allopatric speciation (Allegrucci et al. 2017). This dynamic has been particularly well-documented in orthopterans with low dispersal capabilities, where repeated isolation and reconnection events, driven by glaciation cycles, have resulted in speciation bursts, especially in geographically complex regions such as the Mediterranean and Anatolia (Ortego et al. 2024).

Widespread European genera, like Calliptamus Serville, 1831, a member of the subfamily Gomphocerinae, feature morphologically similar species but display distinct delimitation patterns, including probable cryptic taxa (e.g., within C. siciliae Ramme, 1927 and C. barbarus (Costa, 1836)). However, these patterns within Calliptamus diverge from those observed across other genera in the Acrididae family, despite their similarly broad distribution across Europe.

Another factor impacting geographical diversification patterns is anthropogenic transport and the introduction of non-endemic species. A well-documented case involves Rhacocleis annulata Fieber, 1853, native to Sicily but subsequently transported across large areas of continental Europe through the ornamental plant trade (Barataud 2018). A similar scenario has been proposed for Eupholidoptera smyrnensis (Brunner von Wattenwyl, 1882) in the Aegean islands (Çiplak et al. 2010b) and may have affected numerous other species, both in recent times and historically, altering distribution ranges and complicating phylogeographic reconstructions (Jesse et al. 2011).

The BIN discordance analysis and observed delimitation patterns have revealed several cases of incongruency with currently accepted taxonomy. In order to exclude the misidentification of specimens and the use of obsolete taxonomic units as reasons for these incongruencies as far as possible, MEDOR was initially based exclusively on specimens collected in the field over the last decade (between 2013 and 2022) which were identified by experts directly after collection. The identification was verified after analysing the ML trees using the most recently published dichotomous keys. As taxonomic revisions of orthopterans are continuously being published (e.g., Barat 2007, 2012, 2013; Massa et al. 2023), further updates for many groups were required. Consequently, this dataset conforms to the most up-to-date taxonomy.

The quality control protocol, involving sequence alignment assessment, phylogenetic reconstruction, chromatogram analysis, and sequence identification using BLASTn, has excluded most suspect records potentially representing numts to the greatest extent possible. Additionally, the improvement in PCR primers has contributed to reducing numts. Furthermore, Moulton et al. (2010) noted that the proportion of amplified numts clustering with an incorrect species is minimal compared to the number of orthologs and paralogs that form independent clusters or no clusters at all, suggesting numts occurrence as a minor concern in this analysis.

The influence of Wolbachia bacteria was also examined by analysing each sequence with BLASTn. However, merely detecting the presence of this endosymbiont does not necessarily indicate its impact on the host genome or the occurrence of horizontal gene transfer (HGT) (Hawlitschek et al. 2016). HGT signatures are typically identified and removed from data based on phylogenetic incongruence, patchy distribution, or compositional anomalies (Zhaxybayeva and Doolittle 2011), with this project only adopting the first approach due to the limitations of genome-wide strategies.

Signatures of potential hybridization and introgression have been confirmed in several taxa within this dataset, with notable cases reported among both caeliferans and ensiferans. Hybrids among Chorthippus biguttulus (Linnaeus, 1758) and related species have been generated under laboratory conditions, but hybrid males were found to be behaviorally sterile due to their intermediate courtship songs rejected by all females (Gottsberger and Mayer 2007), suggesting HGT in the wild to be minimal or absent. On the other hand, as reported by Korsunovskaya (2016), all intermediate morphological forms between Platycleis albopunctata (Goeze, 1778) and P. intermedia (Serville, 1838) emitted calling signals that match the temporal pattern of P. albopunctata, indicating that hybrid males may reproduce. DNA barcoding cannot detect genetic admixture, as the mitochondrial COI is exclusively inherited through the maternal lineage. Therefore, the presence of hybrid specimens in the current dataset cannot be ruled out, requiring further evidence for confirmation (Zhaoke et al. 2021).

Incomplete lineage sorting (ILS), the sharing of ancestral haplotypes among related species, is identified as a primary cause of barcoding failure in numerous groups. Rapid and recent speciation, particularly within the Gomphocerinae subfamily of Acrididae, which diverged approximately in the last 6.44 million years (Hawlitschek et al. 2022), and specifically within the genera Chorthippus and Stenobothrus Fischer, 1853 (Vedenina and Mugue 2011), exemplifies this issue. Identical mitochondrial haplotypes and barcodes can be obtained between populations isolated for ~1 million years (Nabholz 2023). Additionally, post-glacial re-colonization has been suggested as the origin for Central European populations of some Acrididae species, potentially leading to recent secondary contacts and narrow hybrid zones (Lunt et al. 1998; Nolen et al. 2020).

The current delimitation analyses underscore the limitations of barcoding in distinguishing extremely recently diverged species, while also highlighting interesting patterns in less recently diverged lineages that may conceal cryptic taxa, due to nearly identical morphological and bioacoustic characteristics among different species.

A prime example is provided by the insular Greek populations of Eupholidoptera Maran, 1953 and Poecilimon Fischer, 1853, particularly species such as E. smyrnensis and P. cretensis Werner, 1903. Despite being classified as a single taxon, populations from different islands or even different regions within the same island (e.g., Crete) exhibit clear genetic divergence. This divergence is influenced by the mutation rate of the marker in question (COI), which dictates the expected number of mutations diagnostic of a species. The mutation rate determines the pace at which neutral divergence accumulates, potentially enlarging the barcoding gap in rapidly evolving species (Allio et al. 2017).

As previously noted, several evolutionary processes, predominantly hybridization and incomplete lineage sorting (ILS), significantly influence the performance of species delimitation based on barcode sequences. Utilizing methods based on diverse theoretical frameworks helps mitigate biases associated with the limitations of each algorithm.

Taxonomically problematic groups, including Platycleidini, Gomphocerinae, and Stenobothrinae, could be more accurately delimited and phylogenetically resolved through the adoption of more rigorous protocols. Schmidt et al. (2024) have shown RADseq to be an effective method for elucidating the phylogeny of the C. biguttulus group and other species within the genera Chorthippus and Pseudochorthippus Defaut, 2012. The species clusters derived from this approach closely correspond with traditional taxonomy, notwithstanding the complexities introduced by their recent evolutionary histories and frequent occurrences of hybridization. Moreover, significant disparities between nuclear and mitochondrial phylogenies were noted, rendering COI alone insufficient for comprehensive phylogenetic analysis and the detection of recent radiative events.

High-throughput genomic techniques promise to effectively identify divergent lineages among both widespread and endemic species. Esquer-Garrigos et al. (2019) uncovered a surprisingly intricate level of biogeographical complexity at a localized scale in populations of Ephippiger diurnus Dufour, 1841 in the Pyrenees. Similarly, the current study has revealed potential divergent lineages within endemic and geographically restricted species, notably among ensiferans, with Antaxius difformis (Brunner von Wattenwyl, 1861) and Thyreonotus bidens (Bolívar, 1887) serving as prime examples.

Future investigations should leverage high-throughput genetic methodologies not only to address existing taxonomic uncertainties but also to explore the biogeographical characteristics of endemic orthopteran species, which may demonstrate greater (rather than less) variation and divergence compared to their more ubiquitous counterparts.