Collecting trips were performed in the mountains of Peloponnese, Central Greece, Evvoia, and the Pindos range in the years 2021, 2022 and 2024, leading to a total of 28 mountain summits visited (Fig. 3). Localities were selected based on the list of mountains from which Parnassiana specimens had previously been collected (Willemse et al. 2018), as well as mountain summits lacking records but exhibiting conditions similar to those of known localities, mainly elevations above 1900 m. Maps were created in QGIS Desktop v.3.28.12 (QGIS.org 2023) using an elevation gradient raster from the SRTM database (Consortium for Spatial Information (CGIAR-CSI) 2004, https://srtm.csi.cgiar.org/srtmdata/), and freely available vector data (IUCN Basedata 2024). Additional outgroup material from taxa of Platycleidini, with a focus on putative phylogenetically close lineages, was collected from Albania, Bulgaria, Kazakhstan, North Macedonia, or downloaded from GenBank (www.ncbi.nlm.nih.gov/genbank).

Details for each specimen used in the present study (coordinates, collection data, locality, GenBank accession numbers, etc.) can be found in File S1. Newly collected specimens are deposited in the Institute of Biodiversity and Ecosystem Research, Bulgarian Academy of Sciences (IBER-BAS) and the Natural History Museum of Crete (NHMC). For the morphological study, a small number of dried specimens from the Naturalis Biodiversity Center’s collection were also examined, including type series and topotypes of several species (P. tymphrestos, P. nigromarginata, P. tymphiensis).

A small set of morphological characters usually used as taxonomically-informative in Platycleidini was examined to complement the molecular analyses. We used these characters to assist with taxa delimitation by differentiating the clades recovered in the molecular analyses and to trace the morphological evolution in a comparative phylogenetic context. Four taxonomically significant characters used in species descriptions are compared here: for males, the last abdominal tergite (dorsal view), the cerci (dorsal view) and the titillator (anterior view); for females, the subgenital plate (ventral view). Morphological terms follow Willemse (1984) and Heller (2006). We consider the term “epiphallus” and “epiphallic sclerites”, that were used in the descriptions of species (Willemse 1973, 1980; Willemse and Willemse 1987; Heller and Willemse 1989), as synonymous to “titillator”.

The Parnassiana specimens collected during the 2021-2024 field trips were preserved in ethanol. For examination they were superficially air-dried, or, for smaller structures, submerged in ethanol, and photographed with a Zeiss Stemi 2000-C Stereo Microscope equipped with a Canon EOS 1300D camera. Genitalia were extracted and, if necessary, treated with a hot 10% KOH solution for a few minutes to remove the soft tissues. The dried specimens studied in Naturalis were placed in humid conditions overnight so the tissues would soften.

Song recordings were made in captivity in the field or in laboratory using a Tascam DR-680MKII (192 kHz, 24 bit) or, in a few occasions, ZOOM H2 (96kHz/24 bit) digital recorders, connected to a Pettersson D500 external microphone, equipped with a custom-made preamplifier. Recorded specimens were usually kept within plastic or paper containers with a volume of 400 ml covered with gauze caps. Microphones were fixed directed to the cages openings at a distance to avoid oversampling. Additional recordings were downloaded from www.xeno-canto.org. Details on specimen collecting sites, recording conditions, and equipment are provided in File S3. Recordings were visualised and analysed in Audacity v. 3.7.7 (https://www.audacityteam.org).

Bioacoustics terminology follows Heller (2006) and Ivković et al. (2017). Calling song—the song produced by an isolated male; syllable—the sound produced by one complete opening and closing movement of the tegmina; microsyllable—a short (impulse-like) syllable; macrosyllable—a long (typical) syllable; echeme—a first-order assemblage of syllables; sequence—syllable/echeme series of different length; phrase—combination of different echemes.

NAD2 was amplified for a total of 91 specimens, COI for 61 specimens, and ITS for 86 specimens. The dataset includes outgroups, all valid Parnassiana species, 27 out of 33 known Parnassiana localities reported in Willemse et al. (2018), and one new locality for the genus (Mt. Koziakas). Details and GenBank accession numbers are given in File S1. Of the three genetic markers, NAD2 has the highest variability, and the nuclear ITS the lowest. Summary statistics for all datasets are shown in Table 1, and single-locus BI trees of the three markers are shown in Files S4–S6. No signs of significant saturation or NUMTs (for the coding fragments) were found in any locus, and ESS values were above 800 for all runs. Single-locus matrices resulted in poorly-supported trees for ITS and COI, although COI mostly recovered the topology of the concatenated matrix, albeit with low support values. The NAD2 trees showed strong support with the exception of the deeper clades, and topology almost identical to the concatenated matrix. As a result of the poorly supported ITS phylogeny and the convergent topology of the strongly-supported mitochondrial markers, the trees based on the concatenated matrices reflect mostly the mitochondrial phylogeny of Parnassiana.

The amplified sequences were combined in a concatenated NAD2–COI–ITS matrix including all available Platycleidini outgroups and Parnassiana individuals with ≥ 2 high-quality genetic markers, which was used for ML and BI analyses and showed high ESS values (above 800 for all runs). The matrix consisted of 69 Parnassiana sequences (ingroup) and 24 outgroup sequences, and had a length of 2863 bp. For divergence time estimation, a concatenated NAD2–COI–ITS matrix which included a single individual from each mountain summit/population was preferred, leading to a total of 45 sequences (31 ingroup and 14 outgroups), and a length of 2862 bp. ESS values for all statistics were high (above 200). For species delimitation, single-locus NAD2 and COI matrices containing only Parnassiana sequences were used, with a size of 67 sequences and length of 993 bp for NAD2, and 47 sequences with a length of 969 bp for COI.

The phylogenetic analyses based on the concatenated matrix (Fig. 2; Files S7, S8) supported the monophyly of the genus Parnassiana. The phylogenetic tree inferred by the BI analysis (Fig. 2) showed higher support and slightly different topology than the ML tree (File S8). The ML tree provides lower support for a monophyletic Parnassiana, and thus accounts for a possible polytomy at the root of the MRCA of Parnassiana and part of Montana. As this suggests a paraphyletic Parnassiana, which contradicts the hypothesis provided by the BI phylogeny, morphology and acoustics, we prefer to base our phylogenetic inferences on the BI tree.

Based on the phylogenetic analyses, three main geographically outlined clades are revealed within Parnassiana: (1) the populations from Peloponnese, except for Mt. Taygetos (P. fusca), henceforth referred to as the Peloponnese clade; (2) the population from Mt. Taygetos, henceforth referred to as the Taygetos clade; and (3) the populations of mainland Greece, henceforth referred to as the Mainland clade (Fig. 3).

The Peloponnese clade is a highly supported monophyletic group consisting of the studied Parnassiana populations from all Peloponnesian summits except for Mt. Taygetos. The topology reflects the relative geographic position of the mountains (Fig. 3): Parnon, the southernmost and most isolated locality of the group, has a basal position; Menalon and Kyllini, in the spatial center of the group’s distribution, form a strongly supported clade with long branches; and the three mountains northwest of the former (Chelmos, Erymathos, Panachaikon) are grouped together.

The position of P. fusca, forming the Taygetos clade, is debatable, as it has low support values (compare Fig. 2 with File S8), but in no analysis it grouped with the other Peloponnesian lineages.

The well supported Mainland clade includes all lineages inhabiting Central Greece, the Pindos range, and Evvoia island (Fig. 3). It consists of four subclades whose relationships were not fully resolved due to low support at the basal nodes: (1) the Akarnanika subclade (P. nigromarginata); (2) the Dirphys-Parnassos subclade (including P. dirphys and P. parnassica); (3) the South Pindos subclade (SPS) – a well-supported lineage that spreads across most mountains of Central Greece (the southern extensions of the Pindos range); and (4) the North-Central Pindos subclade (NCP) – a group inhabiting most mountains of the Central and North Pindos range. The last two further diversify to multiple independent lineages that reveal complex geographical patterns.

The South Pindos subclade (SPS) is mainly distributed in Central Greece and consists of the currently described taxa P. gionica, P. coracis, P. tymphrestos and P. panaetolikon, as well as multiple populations with unknown taxonomic status. P. gionica, known only from Mt. Giona, branches off first. Mount Vardousia, the type locality of P. coracis, and Mt. Oxia, form the next highly supported cluster. Descending from these nodes, there is a monophyletic group including eight mountains. Mount Tymphrestos, the type locality of P. tymphrestos, has its lineage most closely related to the population of Mt. Kaliakouda. The lineage of Mt. Oiti, currently assigned to P. tymphrestos (Willemse 1973), is also related to these, though with lower support. The other branch consists of the lineage from Mt. Panaetolikon, type locality of P. panaetolikon, and a geographically distant but strongly supported cluster of lineages originating from four different mountains (Avgo, Agrafa, Gavrogo, Tzoumerka) (Fig. 3). These four lineages are scattered throughout the summits of the Central Pindos range, very distant geographically from their sister group P. panaetolikon and from all other lineages of the SPS. Despite the distant and disjunct locations of these panaetolikon-like populations (Heller and Willemse 1989; Willemse and Willemse 2008), the genetic distances between them are remarkably low. To sum up, the SPS includes almost all lineages found in the main massifs of Central Greece, as well as a taxon which has spread over a number of summits of the Central Pindos range (Parnassiana sp. 3) (Figs 2, 3).

The North-Central Pindos subclade (NCP) is the least known Parnassiana subclade, as it includes only two lineages affiliated to described species, P. tymphiensis and P. tenuis, and seven lineages that are not ascribed to taxa (compare Willemse et al. 2018). Its monophyly is well-supported, but the deeper relationships between its lineages are not resolved (Fig. 2). A well-outlined cluster within NCP is [(Tymphi + Koziakas) + (Avgo + (Karava + Voutsikaki))]. The other lineages comprising the NCP include five populations with unresolved relationships, yet clearly distinct from each other. These include: P. tenuis from Mt. Tzoumerka (type locality); a population from Mt. Voutsikaki (Parnassiana sp. 5), which is distinct from the individuals from Voutsikaki and Karava presented here as Parnassiana sp. 2; the population from Mt. Smolikas, formerly assigned to P. tymphiensis (Willemse 1980); and the lineages found on Mt. Delidimi and Mt. Helidona. Regarding the spatial distribution of NCP, all lineages but one are found in the Central and North Pindos range, with Mt. Delidimi being the southernmost, and Mt. Smolikas the northernmost locality (Fig. 3). A single lineage, however, inhabits Mt. Helidona in South Pindos, which is in close proximity to the mountains of Panaetolikon and Kaliakouda. Interestingly enough, the nuclear marker (ITS) alone places this population as sister to P. panaetolikon with high support, showing an important discrepancy between the mitochondrial and nuclear markers (see File S6).

The results of two of the calibration schemes – the biogeographical calibration and the NAD2 evolutionary rate for micropterous Orthopterans (Chang et al. 2020) – converged, giving estimates of less than 0.1 million years difference for the deepest nodes and much less for the more recent ones. They are presented on the BEAST chronogram in Fig. 4, while the full tree based on the micropterous molecular rate, including outgroups, can be found in the File S9. Their convergence, in addition to the strong correlation of the dating with climatic events in the area, suggest the resulting dating as the most likely scenario. Posterior probabilities and 95% HPD intervals are shown in File S10.

The estimation based on the evolutionary rate for Pholidopterini (Çiplak et al. 2022) provided younger age estimates than the other methods of up to one million years difference for the deepest Parnassiana nodes. The resulting dates also showed low correlation to climatic events. Therefore, it is treated here as a less likely scenario and will not be further discussed. Time estimates, HPD intervals and posterior probabilities based on the Pholidopterini molecular rate can be found in File S10.

The divergence of Parnassiana from its most closely related lineage of the genus Montana is dated around 4.86/4.76 mya (evolutionary rate/biogeographical calibration) (node 1). The Most Recent Common Ancestor (MRCA) of Parnassiana, marking also the first major split into Peloponnese clade and [Taygetos clade + Mainland clade], is dated around 4.16/4.08 mya (node 2). It is followed closely by the split between the Taygetos and Mainland clades, which was used as the calibration point for the biogeographical calibration, at 3.86/3.79 mya (node 3). The poorly-supported position of P. fusca adds a measure of uncertainty at the timing of these basal nodes, but does not significantly affect the time estimates of the descending nodes, as was found by trials positioning P. fusca as a basal group of Parnassiana (trees not shown).

Descending into the main clades and subclades of Parnassiana, a number of temporal patterns can be observed. The MRCA of the Peloponnese clade (3.06/3.01 mya) and the Mainland clade (3.16/3.1 mya) were dated closely. Inside the Mainland clade, the speciation events that lead to the formation of the four subclades took place in quick succession at the end of the Pliocene, from 3.16/3.1 to 2.86/2.8 mya (nodes 4, 10, 11, 12). A second rapid radiation followed in the Early Pleistocene (2.64/2.59–2.45/2.41 mya), giving rise to the main nodes of NCP. The SPS, on the other hand, diversified later (MRCA was dated 2.24/2.19 mya; node 16), following a more gradual pattern. In total, eight out of a total of 31 cladogenic events, including major clades and subclades, took place during the Late Pliocene–Early Pleistocene (3.16–2.41 mya).

Moreover, time estimates for the majority of the cladogenic events of Parnassiana (12 nodes out of 31) fall within the Early Pleistocene (2.42/2.38 to 0.98/0.96 mya), while only a few (eight nodes, four of which correspond to a single taxon), are of Middle Pleistocene age or younger (0.77/0.75 mya). Lastly, only three lineage clusters (nodes 27–31) show time estimates matching or younger than the Mid-Brunhes Event (0.43 mya). The later data indicate long isolation periods (over one million years) between most Parnassiana populations, while only eight out of the 28 populations studied have had mitochondrial gene exchange with each other in the last 500.000 years.

Sequence-based species delimitation tests, which estimate the gap between intra- and inter-specific distances, were performed on single-locus sequences. ITS was not used due to its poor performance in resolving the phylogenetic relationships of Parnassiana.

A GMYC analysis was performed only on the NAD2 genetic marker, as it was the one with the most sequences and best bootstrap support values. The single-threshold GMYC test results showed a total of 27 ML entities with a confidence interval of 24–31 (File S12). Multiple thresholds tests overestimated the number of species significantly (results not shown).

The ABGD and ASAP methods were tested both on the NAD2 and COI markers independently. Full results are shown in the File S13. All runs of ABGD with X between 0.5–1 on NAD2 showed identical results and suggested 23 putative species, providing the least number of species across all used delimitation methods (Fig. 2). The five partitions with the best ASAP score suggest between 23 and 30 taxa. For the COI marker, all runs in ABGD with X between 0.5–1 showed identical results and suggested 18 species. Excluding the five taxa not presented in the COI dataset, the result corresponds to the ABGD results for NAD2, while for ASAP, the five partitions with the best score suggest between 20 and 25 species excluding the five taxa that were not represented.

The comparison of morphological characters used as diagnostic in species descriptions, including male cerci and titillators (epiphallic sclerites) (Fig. 2), the shape of male 10^th abdominal tergite (Fig. 5), and the female subgenital plate (Fig. 6), revealed that most terminal lineages lacked exclusive morphological synapomorphies, with only major clades recovering common patterns of their titillators and cerci. Τhe taxa belonging to the Peloponnese clade share large titillators with strong and long apical arms shaped in a similar way, and cerci that are stout and wide basally. These characters are mostly absent from the Mainland clade, with the exception of the most basal lineages representing P. nigromarginata (large titillators with strong apical arms and stout cerci) and P. dirphys and P. parnassica (stout and wide basally cerci). P. fusca exhibits intermediate characters, with the 10^th abdominal tergite and cerci being similar to the Peloponnese clade, specifically to P. chelmos, while the smaller titillators lacking dark coloration with shorter apical arms seem to be transient to the Mainland clade.

In-group differences between the members of each major clade are less pronounced. Especially NCP shows a lot of homogeneity in the morphological characters of its members, highlighting cryptic diversity revealed by molecular data. On the other hand, there are a few striking examples of apomorphies both in basal and terminal lineages of the Mainland clade in the shape of cerci (P. nigromarginata, P. tenuis, Parnassiana sp. 2 + sp. 4), titillator (P. nigromarginata, P. tenuis, Parnassiana sp. 6), male 10^th tergite and female subgenital plate (P. tenuis). Pronotum shape and relative length of tegmina also characterize some of the clades and certain lineages (compare Fig. 1). The pronotum disk is smoother and roundish in the basal lineages (Peloponnese and Taygetos clade, Akarnanika and Dirphys-Parnassos subclades). Length of tegmina relative to pronotum is shorter again in the basal lineages except for the Taygetos clade; in the SPS tegmina are longer (in males about the length of pronotum), while in the NCP male tegmina are longer than pronotum and the pronotal disk is distinctly indented from both sides of the medial keel. Dark coloration is typical for the Peloponnese clade.

In the present study we examined the song pattern in almost all known populations and phylogenetic units except for those from Oiti and Helidona. Songs are characterized by phrases consisting of series (echemes) formed by two types of syllables – microsyllable series (consisting of short syllables produced by faster movement of tegmina) and macrosyllable series (consisting of long syllables produced by slower movement of tegmina) (Figs 7, 8). Both microsyllables and macrosyllables are complex as they are produced in pairs by combined lower-amplitude faster and higher-amplitude slower movement. Song characteristics for each putative taxon are provided in File S11.

Phylogenetic lineages within Parnassiana are characterized by the position of the microsyllable series in relation to the macrosyllable series (echemes) (before, after or isolated), and the length of the micro- and macrosyllable series defined mainly by the number of syllables within each echeme. Similarly to the morphological characterization, the basal lineages (Peloponnesian and Taygetos clade, Akarnanika, Dirphys-Parnassos subclades, but also P. gionica) group acoustically by the ancestral positioning of the microsyllable series before the macrosyllable echemes; however, in the lineages from Chelmos, Panachaiko and Erymanthos, as well as possibly in P. parnassica, the position was reversed. Uniquely, in P. nigromarginata, microsyllable series may appear both before and after the macrosyllable echemes. In the studied terminal lineages of the SPS and NCP subclades, apart from P. gionica, microsyllables are produced only after the macrosyllable echemes. In agreement with Heller’s (2006) statement, both SPS and NCP show significant variation in the duration (number of syllables) of the macrosyllable echemes, that may occur within the same individual.

The present study is a first attempt to reconstruct the phylogenetic relationships of the genus Parnassiana based on a combined set of one nuclear and two mitochondrial DNA fragments.

ITS has been previously used in Orthoptera as a nuclear marker for species and genus-level phylogenetic reconstructions (Kaya et al. 2013; Çiplak et al. 2015, 2020; Borissov et al. 2021, 2023; Kociński et al. 2022), instead of the more commonly used Histone 3, Wingless, 28S and 18S rDNA, which have much lower mutation rates and are usually preferred for higher-level phylogenies (Simon et al. 2006; Song et al. 2015, 2018; Mugleston et al. 2018). In Parnassiana, however, extensive deletions and insertions in ITS1 and ITS2, combined with a low number of parsimony informative sites, have led to low resolution and large polytomies from which few insights can be drawn, and thus the influence of the ITS matrix on the concatenated tree was limited (File S6). Similarly, low resolution of ITS for species-level phylogenies was observed in other Tettigoniinae from the tribe Pholidopterini (Çiplak et al. 2020, 2022), contrary to its high phylogenetic performance in Phaneropterinae as shown for Barbitistini (Chobanov et al. 2017; Borissov et al. 2023).

The mitochondrial marker NAD2 provided the highest resolution and statistical support among the markers (File S4), and had the most significant influence on the topology of the concatenated gene tree. NΑD2 is known to be excellent for revealing species-level phylogenies in Orthoptera, as it has a high proportion of variable and informative sites and provides well-resolved trees (Simon et al. 1994, 2006; Çiplak et al. 2015, 2020; Borissov et al. 2021; Kociński et al. 2022; Willemse et al. 2023). COI, a marker very common in comparative and species-delimitation studies, provided similar topology to NAD2, albeit with fewer parsimony-informative and variable sites and lower support values (File S5). Such results indicating that NAD2 has more parsimony-informative and variable sites than COI and higher branch support are in line with previous observations in Orthoptera genera, such as Poecilimon and Isophya (Chobanov et al. 2017; Borissov and Chobanov 2020), and in Odonata (Cheng et al. 2018).

Comparison of the mitochondrial and nuclear phylogenies can reveal signs of mitonuclear discordance. Due to the poor resolution of the ITS tree, however, the concatenated gene tree largely reflects the mitochondrial evolutionary history of Parnassiana, and we rely on signs from phenological characters, mostly those subjected to sexual selection, i.e., genitalia and acoustic communication, to detect such discordances.

The phylogenetic tree presented here is mostly congruent with the current systematics of the genus (Massa and Fontana 2011), as it supports the monophyly of Parnassiana and all thirteen nominal species. Furthermore, the phylogenetic analysis and species delimitation tests support the observation by Willemse and Willemse (2008) that certain populations found in isolated massifs constitute unnamed species.

Signs of discrepancies between the presented gene tree and species delineations based on morphology are observed in several cases. Τhere is at least one clear example of mitochondrial introgression in our phylogeny. Within the SPS, Parnassiana sp. 1 groups with P. tymphrestos but shows clear morphological and acoustic characteristics of the P. panaetolikon + Parnassiana sp. 3 group. In the Peloponnese clade, the lineage from Mt. Kyllini (P. chelmos unicolor), formerly grouped under the same taxon as the Panachaikon population (Willemse 1973, 1980), is sister to P. menalon instead of the other populations currently placed within P. chelmos (Willemse 1973, 1975, 1980). The molecular phylogeny is further supported by the song pattern of the Kyllini population (Fig. 7C), which clearly groups with P. parnon and P. menalon based on the long phrases starting with microsyllable series, while in other members of P. chelmos the microsyllable series are produced at the end of the phrases (Fig. 7D–F). Yet, the Kyllini lineage shows clear morphological similarity with P. chelmos (Figs 1, 2, 5, 6). There are two possible scenarios to explain this pattern. One scenario is that the Kyllini population kept mitochondrial DNA from an extinct sister species to P. menalon that once existed on Kyllini, whose gene pool was flooded by immigrating males from the surrounding mountains hosting P. chelmos, a process called ghost introgression (Shen et al. 2025). This scenario requires a convergent evolution of the song pattern in Kyllini males that may be driven by sexual selection of the females of the local species. Another explanation of the observed phenomenon may be that the morphological characters of P. chelmos are plesiomorphies shared by the entire clade, whereas the unique morphology of P. menalon presents an apomorphy or atavistic reappearance (a scenario that is partly supported by the ITS phylogeny – File S6). Another intriguing case which could be related to ghost introgression is Parnassiana sp. 6 (Mt. Helidona). Based on the mitochondrial phylogeny, the population belongs to the NCP, but geographically it belongs in Central Greece and shows morphological similarity with the P. panaetolikon and Parnassiana sp. 1 from its neighboring mountains Panaetolikon and Kaliakouda. Moreover, the ITS data clearly group it with P. panaetolikon.

Such discrepancies suggest that the mechanisms that determine the evolution of Parnassiana are complex, and a simple allopatric speciation model may not be sufficient to interpret them, a pattern that is common among sky island inhabitants (Recuero et al. 2014; Knowles and Massatti 2017; Ortego and Knowles 2022).

Based on such an interesting phylogenetic pattern, formerly discussed phenological peculiarities in the genus may be discussed in an evolutionary context. Former studies supported the view that speciation processes in Parnassiana are significantly influenced by sexual selection on male genitalia (titillators, last tergite, cerci) and less on the evolution of calling songs, which are similar between most described taxa and variable within species, therefore having poor systematic value (Heller 2006). Indeed, morphology characterizes the main phylogenetic lineages, with the Peloponnese and Taygetos clades, and the Akarnanika and Dirphys-Parnassos subclades showing clear morphological differentiation. Additionally, a few outliers of the mostly uniform morphology with two types of titillators (short- and long-armed) observed at the terminal nodes of the SPS and NCP, can be observed in P. tenuis (various characters) and Parnassiana sp. 4 (in cerci), both sympatric with Parnassiana sp. 3.

Song characteristics also support the main lineages, though their pattern is different from the one observed in the morphology, suggesting a distinct evolution from it. Here, main song-type groups outline the following groupings: 1. Parnon+Menalo+Kyllini+Taygetos (long to medium macrosyllable echemes with microsyllable series produced isolated or before the macrosyllable series); 2. Chelmos+Panachaiko+Erymanthos (short to medium-length macrosyllable echemes with long microsyllable series produced after the macrosyllable series); 3. Akarnanika (unique pattern with microsyllable series produced both before and after the macrosyllable series); 4. Dirphys (short macrosyllable echemes with microsyllable series produced before the macrosyllable series containing long macrosyllables); 5. Parnassos and Giona (two types of songs with very long nocturnal macrosyllable echemes but differing in the production of microsyllables); and 6. the rest of the SPS together with NCP (variable short to medium-length macrosyllable echemes with short microsyllable series produced after the macrosyllable series). It is worth noting the resemblance of song types of distantly related lineages like P. chelmos, the SPS (except P. gionica), and the NCP, that obviously evolved independently superficially similar song pattern apomorphies.

Significant intraindividual variation was observed within the SPS and NCP subclades depending on the ambient temperature, with specimens tending to produce denser and shorter (composed of less syllables) macrosyllable echemes with increased temperature (examples from own recordings with known details shown in Fig. 8J1 and J3, M1 and M2). At the same time, continuously singing animals tend to change the temporal song pattern with time, which may reflect change in the body temperature (warming up) at similar ambient temperature (Fig. 8C1 and C2, E1 and E2). Yet, these variations do not reach the scope of differentiation in song patterns between sympatric taxa, i.e., Parnassiana sp. 3 (Fig. 8G) occurring together with Parnassiana sp. 4 (Fig. 8H) and P. tenuis (Fig. 8K).

Regarding the sympatric taxa, Heller and Willemse (1989) suggested that the songs of two taxa found in the same location (P. tenuis and an unnamed species in Tzoumerka) show very limited differences in song pattern compared to their morphological distinction. However, according to our observations, the P. tenuis song considerably differs from the syntopic Parnassiana sp. 3, having much shorter macrosyllable echemes (compare Figs 8G and K). Similar song differentiation was observed in the syntopic Parnassiana sp. 3 and sp. 4 on Mt. Avgo (Figs 8G and H).

Much more interesting are the observed examples of different echeme length depending on the time of the day. Although there are no distinct diurnal and nocturnal song patterns in Parnassiana (contrary to some related groups like Montana; see Ivković et al. 2017), the number of syllables within the macrosyllable echemes and thus their length may differ when songs are produced in day or night (see Fig. 7J, K for examples of diurnal and nocturnal songs produced at similar temperature). In P. parnassica and P. gionica, nocturnal songs (Fig. 7J3, K3) contain many more macrosyllables within the series than diurnal songs (Fig. 7J1, K1) at similar temperature. In those cases, temperature seems not to be a major factor and thus the case may concern a primitive stage of distinction between diurnal and nocturnal song patterns known in some other related genera (e.g. Ivković et al. 2017; Barataud 2025).

The above discussed phenotypic evolution of morphological and acoustic traits of Parnassiana supports the complex phylogenetic pattern and proves that speciation processes in the group have been led by sexual selection that accelerated diversification in sympatric taxa. Both the evolution of genitalia and behavior involved in reproduction show strong phylogenetic signals but followed distinct evolutionary pathways.

Parnassiana originated in the southern areas of the Balkan Peninsula ca. 4.86/4.76 mya in the Middle Pliocene from a common ancestor shared with a yet poorly outlined infragroup of Montana (Fig. 4). Shortly after, the Corinthian Rift developed, causing the isolation of Peloponnese from mainland Greece 4–3.5 mya (Wicker et al. 2024), and the first major split of Parnassiana into Peloponnesian clade + [Taygetos + Mainland clades] (4.16–4.08 mya). Around the same time (4–2 mya) the uplifting of Mt. Taygetos occurred, which, together with rifting that formed deep and occasionally flooded valleys, led to the massif’s isolation from the rest of the Peloponnese (Fountoulis 2014; Kleman et al. 2016) and to the permanent isolation of P. fusca both from the highly-supported Peloponnesian clade and from the Mainland clade. These two processes resulted in the geographic split of the ancestral lineages of Parnassiana and the establishment of its three main clades (Peloponnesian, Taygetos and Mainland) via vicariance. This timeline of divergence 4 to 3 mya is shared by many different organisms in the area, such as the land snail Codringtonia (Kotsakiozi et al. 2012), the mammal Talpa stankovici (Colangelo et al. 2010), lizards (Psonis et al. 2018), the scorpion genus Euscorpius (Parmakelis et al. 2013), and a number of Dolichopoda species (Allegrucci et al. 2009, 2021).

In the Late Pliocene (3.6–2.58 mya; Jiménez-Moreno et al. 2013), constant climate cooling and drying transformed the savannas and subtropical humid forests of the area to dominating sclerophyllous and steppe vegetation with the establishment of drier seasonal Mediterranean climate 3.2 to 2.8 mya (Suc 1984). This process has been shown to cause extensive diversification across the Mediterranean in Orthoptera (Çiplak et al. 2015; Allegrucci et al. 2021; Uluar et al. 2023) and other organisms (Kotsakiozi et al. 2012; Fiz-Palacios and Valcárcel 2013; Poulakakis et al. 2015; Kougioumoutzis et al. 2021; Koutroumpa et al. 2021). This cooler and drier climate would have promoted wider establishment of the cold-adapted Parnassiana in the high altitude mountain ranges of the southern Balkan Peninsula. At the same time, long-distance dispersal and contact between populations was prevented due to the fragmentation of the most isolated mountain massifs by a network of lowlands with unsuitable climate. As a result, during this period, migration and isolation events may have promoted secondary diversification within the Peloponnese clade and the establishment of the four major subclades within the Mainland clade. In fact, eight out of a total of 31 cladogenic events, including major clades and subclades, took place during the Late Pliocene–Early Pleistocene (3.16–2.41 mya), marking this time period as a major time stamp in the evolution of Parnassiana, during which the major clades and subclades of the genus were established and the NCP showed a burst of early diversification. In the Pleistocene, climate continued to change and climate cycles that marked the evolutionary pattern of the European biota were established (Hewitt 1996, 2004; Taberlet 1998; Wallis et al. 2016). During the Early Pleistocene (2.58–0.77 mya), short glacial periods of low intensity alternated with interglacials at 41 kyr cycles (Willeit et al. 2019; Shackleton et al. 2023). The majority of extant Parnassiana lineages evolved during this period (Fig. 4), possibly as a result of short stepping stone dispersals during the mild glacial periods followed by isolation periods preventing gene flow between the populations due to the warm interglacials. By the Mid-Pleistocene Transition 1.25–0.7 mya, which marks the end of this stage (Willeit et al. 2019; Shackleton et al. 2023), all putative Parnassiana species were established.

The Mid-Pleistocene Transition marked a significant change in Earth’s climate from “mild” to “full” glacial periods (Willeit et al. 2019; Shackleton et al. 2023). Glacial-interglacial cycles became longer, from 41.000-year to 100.000-year periods, with shorter interglacials, and longer and more intense glacial periods (Barth et al. 2018; Shackleton et al. 2023). Nevertheless, the intense glacial-interglacial cycles of the Middle and Late Pleistocene appear to have had limited effects on the speciation of Parnassiana. Only a small number of lineages diversified after the Mid-Pleistocene Transition (e.g., P. coracis, P. chelmos), and these splits produce only intraspecific divergence. Only three taxa show signs of recent long-distance dispersal (either genetic and/or geographic – Parnassiana sp. 3, P. tymphiensis and P. tenuis), taking advantage of the intense glacial periods to expand their ranges using the Pindos mountains as a corridor.

The evolutionary history of Parnassiana is reflected in its geographic patterns. Massifs isolated by larger distance and low-altitude valleys, such as Taygetos, Dirphys, and Akarnanika, host well-defined species with long evolutionary histories, reflected in their morphological and song apomorphies. On the other hand, summits that are parts of large massifs, such as the mountains of northern Peloponnese, the Central Greece mountains and the Pindos range, are characterized by taxa with more uniform phenotype and complex phylogenetic relationships, with some exceptions that may represent splits of currently isolated or extinct lineages. Such patterns have been previously observed for inhabitants of sky islands such as grasshoppers, spiders and birds (Masta 2000; Robin et al. 2010; He and Jiang 2014; Ortego and Knowles 2022). Secondary dispersals of Parnassiana are characteristic for the Pindos range and are facilitated by the extensive high-altitude range and thus good connectivity and cooler climate. Populations here tend to have wider distributions and higher densities than populations inhabiting isolated small mountain summits (like P. parnon, P. menalon, P. nigromarginata, P. dirphys; own observation). Such dispersals and secondary contacts resulted in the complex phylogenetic pattern of the crown clades: populations from distant mountains grouping in the tree (P. panaetolikon and Parnassiana sp. 3); neighboring mountains hosting distant lineages (Helidona and Kaliakouda); and distinct lineages with syntopic occurrence (Tzoumerka, Avgo, and Voutsikaki).

The clear correlations between cladogenic events, geographical distribution and geo-climatic events, suggest that allopatric speciation was the predominant, though not exclusive, mechanism shaping the evolution and speciation of Parnassiana. This is a very common pattern observed in sky islands (He and Jiang 2014; Martinez-Sañudo et al. 2022), and is primarily driven by niche conservatism, “the tendency of lineages to maintain their ancestral ecological niche” (Wiens 2004). Niche conservatism seems to be especially important in Parnassiana (Stefanidis et al. 2024, 2025), as is the case for many montane species (He and Jiang 2014). Close relatives of Parnassiana, such as Montana and Modestana, are typical of cool climates with most species occurring in mountains or cool deserts (Massa and Fontana 2011). It is therefore expected that the ancient taxon that gave rise to Parnassiana would also share these climatic preferences, and inherit them to its descendants.

Parnassiana populations diversified during the Late Pliocene and Early Pleistocene. The old divergence times indicate that populations, with a few exceptions (Parnassiana sp. 3, P. chelmos, P tymphiensis, P. tenuis), survived locally in the sky islands that acted as interglacial refugia (see Berger et al. 2010 and Çiplak et al. 2015 for similar scenarios), without successfully colonizing nearby mountains, through several glacial cycles. A remarkably similar evolutionary pattern and history was observed in Anterastes, a bush-cricket genus with sky-island distribution in Anatolia (Çiplak et al. 2015; Uluar et al. 2023), highlighting the importance of Pliocene–Early Pleistocene events in the Orthoptera inhabiting the sky islands of the Eastern Mediterranean.

Species delimitation analyses based on the NAD2 and COI molecular markers suggest a range of 18–41 taxa in Parnassiana (Fig. 2; Files S12, S13). ABGD based on the complete sample of NAD2 provided the most conservative number of species with 23 putative taxa, which is congruent with the lowest of the five best estimates of the NAD2-ASAP tests. The NAD2-based GMYC analysis, however, suggested a larger number of taxa (27). GMYC is known for overestimating the number of species, especially when dealing with a small number of individuals per species and strong within-species divergence, which is frequent in groups with geographically isolated populations like Parnassiana (Luo et al. 2018; Hofmann et al. 2019). Additionally, species delimitation analyses based on single-locus DNA do not account for gene tree vs species tree discordance, incomplete lineage sorting, introgression, or gene flow (Luo et al. 2018; Puillandre et al. 2021), processes which appear to have affected the evolution of Parnassiana. Taking into account these inherent biases and the studied phenotype characteristics, we expect that the true number of taxa falls within the lowest estimates of the delimitation tests.

A few notable discrepancies between species delimitation tests and phenotypic characters are discussed below. The ABGD test on the NAD2 marker suggested the lumping of the population of Mt. Kaliakouda with P. tymphrestos from Mt. Tymphrestos, even though the titillators of the two populations bear remarkable differences, and GMYC and ASAP results indicate that the two are distinct species (Fig. 2). As a result, we suggest that Mt. Kaliakouda hosts a yet undescribed species of Parnassiana. On the opposite end, even though all methods recover the specimens from the mountains hosting P. chelmos (excluding Mt. Kyllini) as two or even three distinct entities (Fig. 2), they can be regarded as subspecies based on their conserved phenotype, as has already been proposed (Willemse 1973; Willemse 1980). Similarly, despite GMYC and ASAP suggesting that the lineage pairs inhabiting the neighboring mountains of Vardousia and Oxia, as well as Karava and Voutsikaki, are distinct entities, we consider them as infraspecific based on the ABGD results, the recent divergence times (younger than 0.6 mya), and the lack of discernible phenotypic differentiation. Therefore, we accept the placement of the Parnassiana population of Mt. Oxia in P. coracis, as suggested by Stefanidis et al. (2025), and suggest that the lineages of Karava and Voutsikaki constitute a single undescribed species.

Overall, based on the present results including phylogenetic relationships, species delimitation tests, and a combination of phenotypic characters, we suggest the occurrence of 23 putative species of Parnassiana in Greece.