Specimens were collected at eight sampling sites in Medvednica Nature Park between June 2021 and March 2022 (Suppl. material 1, Fig. 1), with the permission of the Ministry of Economy and Sustainable Development of the Republic of Croatia (UP/I612-07/21-48/73). Adults were collected using an entomological net, while larvae were collected by hand under stones and branches. All collected individuals were preserved in 96% ethanol and stored at +4°C.

Figure 1.  

Sampling sites on the map of Medvednica Nature Park.

Morphological determination of adults was conducted using a ZEISS SteREO Discovery V20 stereomicroscope along with various identification keys (Illies 1955, Kaćanski and Zwick 1970, Sivec and Stark 2002, Zwick 2004, Murányi 2011). The characters used to distinguish species included the patterns of colouration of the head and pronotum, as well as the genital apparatus in both males and females. Larval identification was conducted through sequence comparison with entries available in databases.

Total genomic DNA was isolated from the tissues of one to three legs of adults and larvae, depending on the specimen's size, using the GenElute Mammalian Genomic DNA Miniprep Kit (Sigma-Aldrich, Germany), following the manufacturer's instructions. To amplify the DNA barcode region, we utilised the universal primer set: LCO-1490 and HCO-2198 (Folmer et al. 1994). PCR reactions and subsequent purifications were carried out following the protocols outlined in Hlebec et al. (2022). Sequencing was performed by Macrogen Inc. (Amsterdam, Netherlands).

Chromatograms were analysed and edited in Geneious Prime 2022.1 (Biomatters, Auckland, New Zealand). Misread nucleotides were manually corrected. The presence of stop codons in all sequences was checked using Mesquite (Maddison and Maddison 2021). The BOLD Identification Engine (Ratnasingham and Hebert 2007) and NCBI BLAST search were used to confirm the authenticity of the amplified product.

For phylogenetic and phylogeographic analyses, we selected species with a broader distribution and a larger number of COI sequences available in the BOLD and GenBank databases. The selected species were Brachyptera seticornis, Isoperla grammatica, Leuctra braueri and Leuctra prima. Downloaded sequences were aligned using the ClustalW multiple alignment tool (Larkin et al. 2007) and poor-quality sequences, sequences shorter than 500 bp or those containing stop codons, were excluded from further analyses. Subsequently, sequences were collapsed into unique haplotypes using the online service FaBox v.1.61 (Villesen 2007). Perla pallida Guérin-Méneville, 1838 (CROPL090-19), Isoperla grammatica (CROPL154-21) and Nemoura avicularis Morton, 1894 (CROPL015-21) were used as outgroups. Phylogenetic analyses were performed using the Maximum Likelihood (ML) approach in RaxML-HPC v.8 on XSEDE (Stamatakis 2014) through the CIPRES online service (Miller et al. 2010). Rapid bootstrap analysis was performed using the GTRGAMMA model with 1000 iterations. Nodes that received bootstrap values (RBS) higher than 70% were considered as supported. The resulting phylogenetic trees were visualised and edited using FigTree v.1.4.4. (Rambaut 2018) and iTOL (Letunic and Bork 2021). To investigate concordance between morphological and molecular species clustering, we applied different species delimitation methods: ASAP (Puillandre et al. 2020) using both p-distance and Kimura-2-Parameter (K2P) distance options, ABGD (Puillandre et al. 2011), bPTP (Zhang et al. 2013) and the BIN (Barcode Index Number) system. Additionally, the ranges of intraspecific uncorrected pairwise distances (p-distances, hereafter referred to as "genetic distance") were calculated for each selected species using the programme MEGA 11 (Tamura et al. 2021). TCS phylogenetic networks (Clement et al. 2002) were created using the statistical parsimony method in the PopART programme (Leigh and Bryant 2015). Geographical distributions of haplotypes were made using the Cartopy package v.0.19 in Python v.3.8.

The ASAP analysis, using both p-distance and K2P distance options, grouped haplotypes into three distinct groups. Conversely, the ABGD analysis, utilising the K2P distance option, grouped all haplotypes into a single group, as did the bPTP analysis. Notably, all analysed sequences were assigned to the same BIN (AAY5851) (Fig. 2). Specimen CROPL392-22 from this study was nested within Croatian specimens. The three groups identified by ASAP were also evident in the phylogenetic network (Fig. 2). Samples from Zrinska gora Mountain grouped with those from south-eastern Europe, Bulgaria and Serbia, while Medvednica Mountain samples grouped with those from central and western Europe (G3). The intraspecific genetic distances for Brachyptera seticornis ranged from 0.15% to 2.28% (Suppl. material 1).

Figure 2.  

Phylogeographic analysis of Brachyptera seticornis, based on the COI gene fragment, 658 bp in length. A ML phylogenetic tree. The terminals correspond to the BOLD IDs. Individual collected in the present study is marked with an asterisk. Bars in different colours represent the results of species delimitation methods; B Phylogenetic network obtained by the statistical parsimony (TCS network). Colour coding corresponds to the groups in the species tree; C Geographical distribution of analysed haplotypes. Colour coding corresponds to the groups in the species tree and the TCS network.

The ASAP analysis, using the p-distance option, grouped haplotypes into four distinct groups (with lineage CROPL389-22 nested within G1). In contrast, the K2P distance option yielded three groups. ABGD analysis, also using the K2P distance option, grouped haplotypes into two groups, while bPTP analysis grouped all haplotypes into a single group. Sequences were assigned to two BINs (AAJ2415, AES9453), with specimens from this study belonging to a unique BIN. The species tree showed two clades, one nested within the other. Two groups identified by ABGD analysis were evident in the phylogenetic network (Fig. 3). Samples from Medvednica Mountain formed a distinct group (G1), while another group (G2) included samples from central Europe (Germany, Austria, Italy and Switzerland) (Fig. 3). The intraspecific genetic distances for Leuctra braueri ranged from 0% to 2.46% (Suppl. material 1).

Figure 3.  

Phylogeographic analysis of Leuctra braueri, based on the COI gene fragment, 658 bp in length. A ML phylogenetic tree. The terminals correspond to the BOLD IDs. Individuals collected in present study are marked with asterisks. Bars in different colours represent the results of species delimitation methods; B Phylogenetic network obtained by the statistical parsimony (TCS network). Colour coding corresponds to the groups in the species tree; C Geographical distribution of analysed haplotypes. Colour coding corresponds to the groups in the species tree and the TCS network.

Morphological examination of the Leuctra specimen (CROPL384-22) from Medvednica Mountain was inconclusive. It showed similarities to several species not represented in the BOLD and GenBank databases: Leuctra prima, Leuctra signifera, Leuctra dalmoni and Leuctra carpathica. Therefore, more detailed phylogenetic analyses, including sequences of those closely-related species, were omitted. COI sequence of collected specimen (CROPL384-22) showed the highest similarity to an undescribed Leuctra species (labelled as Leuctra sp. ZB, from Žumberačko gorje Mountain in Croatia, CROPL248-21) (Hlebec et al. 2022). Hence, we adopted the same designation here. Both of these lineages exhibited the closest morphological similarity to Leuctra prima.

The species tree was characterised by two singletons and two highly-supported clades labelled as G1–G4 (Fig. 4), corresponding to two morphotypes: Leuctra sp. ZB and Leuctra prima. All species delimitation methods were congruent in delimiting four groups, except for bPTP, which clustered all Leuctra prima sequences and separated Leuctra sp. ZB (CROPL384-22 and CROPL248-21) into two distinct groups. Additionally, a unique BIN (AES3404) was assigned to specimen labelled Leuctra sp. ZB. (CROPL384-22) collected in this study.

Figure 4.  

Phylogeographic analysis of Leuctra prima and Leuctra sp. ZB, based on the COI gene fragment, 658 bp in length. A – ML phylogenetic tree. The terminals correspond to the BOLD IDs. Individuals collected in present study are marked with an asterisk and Leuctra sp. ZB specimens with triangles. Bars in different colours represent the results of species delimitation methods. B – Phylogenetic network obtained by the statistical parsimony (TCS network). Colour coding corresponds to the groups in the species tree. C – Geographical distribution of analysed haplotypes. Colour coding corresponds to the groups in the species tree and the TCS network.

The four groups identified by ASAP, ABGD and BIN assignment were confirmed in the phylogenetic network (Fig. 4). Interspecific genetic distances between Leuctra prima (G3 and G4) and Leuctra sp. ZB (G1 and G2) ranged from 12.0 to 13.5%, while intraspecific genetic distances for Leuctra prima (G3 and G4) and Leuctra sp. ZB (G1 and G2) reached a maximum of 3.65% (Suppl. material 1).

ASAP analysis, using both the p-distance and K2P distance options, grouped haplotypes into six groups. ABGD analysis, using the K2P distance option, resulted in five groups, while bPTP analysis grouped haplotypes into three groups. Seven BINs were assigned: AAY9655, AEH6396, AEG4373, AEC9627, ACJ0709, AAK4351 and AER8749, with all specimens from this study falling into one distinct BIN. The species tree was characterised by four clades labelled as G1–G3 and G5–G6 (Fig. 5) and group G4 which indicates several singletons. The six groups identified by the ASAP analysis were also evident in the phylogenetic network (Fig. 5). Intraspecific genetic distances for Isoperla grammatica ranged up to 11.7% (Suppl. material 1).

Figure 5.  

Phylogeographic analysis of Isoperla grammatica, based on the COI gene fragment, 658 bp in length. A – ML phylogenetic tree. The terminals correspond to the BOLD IDs. Individuals collected in present study are marked with asterisks. Bars in different colours represent the results of species delimitation methods. B – Phylogenetic network obtained by the statistical parsimony (TCS network). Colour coding corresponds to the groups in the species tree. C – Geographical distribution of analysed haplotypes. Colour coding corresponds to the groups in the species tree and the TCS network.

Brachyptera seticornis, typically absent in Italy, inhabits southern and central Europe, favouring cold streams at altitudes from 300 to 2900 m, including the Alpine Region (Schmidt-Kloiber and Hering 2019).

Species delineation methods revealed up to three groups within B. seticornis, with intraspecific genetic distance up to 2.43%. The haplotypes were separated by a low number of mutations (Fig. 2). One specimen from Medvednica Mt. (CROPL392-22) aligns with samples from Germany and Switzerland (G3), while others from Croatia cluster with those from Bulgaria and Serbia (G1). The population (G2) between central (G3) and south-eastern Europe (G1) comprises the remaining Croatian individuals from Medvednica Mt. and Zrinska gora Mt. Based on the results of the phylogeographic analyses, we assume that the dispersal of the species was from southeast to central Europe (Fig. 2).

The genus Leuctra Stephens, 1836 is a dominant member of the Leuctridae family, characterised by small, dark-coloured species. Resolving phylogenetic relationships within Leuctra remains challenging, with overlapping intraspecific and interspecific genetic distances (Vitecek et al. 2017). To enhance our understanding, we analysed three species from the region: Leuctra braueri, Leuctra prima and an undescribed species, labelled as Leuctra sp. ZB.

Both L. braueri and L. prima are found primarily in mountainous regions of central Europe, particularly the Carpathians, at altitudes above 200 m. L. braueri extends into the Dinaric Western Balkans, while L. prima inhabits the entire Balkans. These species favour cold streams (6 – 10°C) with neutral pH and either moderate to high (L. prima) or low (L. braueri) flow (Schmidt-Kloiber and Hering 2019).

Our discovery of Leuctra braueri marks its first recorded presence in Croatia. It shows genetic divergence, with central European samples (G2) distinct from Croatian specimens (G1). Fewer mutational steps and low intraspecific genetic distances (1.84–2.46%) could suggest a relatively recent intraspecific split. The phylogenetic network hints at Croatian samples potentially being the diversification centre of the central-European lineage, suggesting a migration route from east or south Europe to the west. Furthermore, the question is, would that intraspecific structuring exist if samples from Slovenia and Austria had been included in the analysis? This question emphasises once again the consequences of undersampling in barcoding and phylogeographic studies.

We confirmed the genetic distinctiveness of the undescribed species, Leuctra sp. ZB, from its closest relative, L. prima (interspecific genetic distances were 12.02–13.55%). This exceeds the average interspecific genetic distances in Plecoptera, which is 11.56% (Zhou et al. 2010). All species delimitation methods separated two Leuctra sp. ZB samples (one from Medvednica Mt. and one from Žumberačko gorje Mt.) into distinct groups. In addition to the morphological similarities with L. prima, the specimen from Medvednica Mt. showed a resemblance to L. dalmoni, L. carpathica and L. signifera, which are distributed in central Europe and the Alps (Vinçon and Múranyi 2007). However, their emergence periods do not align with our specimen's late summer collection, as L. carpathica and L. signifera emerge during autumn and L. dalmoni emerges during winter and spring. Furthermore, the distribution ranges of L. carpathica and L. signifera exclude the Balkans and the Dinaric Region, as they are only recorded in the Alps (L. carpathica and L. signifera) and the Carpathians (L. carpathica) (Schmidt-Kloiber and Hering 2019). Regarding the other DNA barcoded Leuctra sp. ZB specimen from Žumberačko gorje Mountain, while the altitude range matches (between 200 and 800 m), the emergence periods do not, as this specimen was found during autumn. Although a longer emergence period is possible, the significant genetic differences between the samples should be considered. Collected specimens likely represent one or two endemic species, with the type localities in north-western Croatia. Detailed morphological examination, including analysis of additional genetic markers, will be necessary to distinguish and describe these species and determine their phylogenetic relationships. The challenge in determining relatedness within the genus Leuctra lies in the lack of genetic data for described species (e.g. Leuctra dalmoni, Leuctra carpathica and Leuctra signifera). It is regrettable that even the most recent study (Reding et al. 2023) did not incorporate molecular data in the description of Leuctra papukensis Reding, Vinçon & Graf, 2023.

Isoperla grammatica is widely distributed across Europe, inhabiting altitudes from lowland to subalpine regions (0–2400 m) with no specific preference for water temperature or pH conditions (Schmidt-Kloiber and Hering 2019).

Phylogenetic reconstruction of Isoperla grammatica revealed complex structuring. Delimited groups had distinct geographical distributions: Croatia (G6, G2), Switzerland (G4), Italy (G1), Austria (G3) and a group (G5) that includes samples from the Alpine Region and central Europe (Slovenia, Italy, Austria and Germany). Croatian individuals formed two separate lineages corresponding to different altitudes: a mountain lineage in the Medvednica Nature Park (G2, altitudes between 300 and 1000 m) and a lowland lineage (G6, altitudes below 200 m). Intraspecific genetic distances ranged from 9.02–11.7%, indicating significant genetic diversity. Expanding sampling to broader scales may increase intraspecific variation, while interspecific divergence impact is likely to decrease due to allopatric speciation dominance (Bergsten et al. 2012). Alpine Region samples were genetically closer to the lowland-Croatian lineage than to Medvednica Nature Park samples, suggesting exceptional isolation of the latter. Genetic distances between Medvednica Mt. population and others ranged from 7.93 to 10.4%, exceeding the maximum intraspecific genetic distances within the genus Isoperla (7.82%) (Hlebec et al. 2021), possibly indicating an undescribed species.

The widespread distribution of I. grammatica in Europe, along with its complex genetic differentiation, suggests the presence of cryptic diversity. Cryptic species lack morphological differences, but are often genetically distinct, forming highly-supported monophyletic lineages (Bickford et al. 2007). Following phylogenetic analyses, individuals are often re-examined to identify distinguishing morphological traits (Pauls et al. 2010). To gain a clearer understanding of the phylogenetic relationships and the delimitation of potential cryptic species within the I. grammatica species complex, a detailed morphological examination, including different life stages and eggs, supported by a multi-gene approach is required. The very recent genomic sequencing of a specimen from the River Test, Great Bridge, Hampshire (McSwan et al. 2023) will be used as a reference for future phylogenetic studies using new-generation techniques.

Stoneflies, known for their poor flying abilities (Stewart 2009), tend to remain in specific areas, resulting in high rates of endemism, mutation accumulation and allopatric speciation. Phylogenetic and phylogeographic analyses of selected species recorded in the studied area (Brachyptera seticornis, Isoperla grammatica, Leuctra braueri) revealed high genetic diversity, except for Leuctra prima.

Phylogenetic reconstructions suggest that central European specimens of Brachyptera seticornis likely originated from eastern Europe, with a population on Medvednica Mt. bridging central and southeast Europe. A similar phylogenetic tree topology was observed in Leuctra braueri, with specimens from Medvednica Mt. representing the first branch-off in species tree.

Individuals of Leuctra braueri, Isoperla grammatica and, potentially, the new species Leuctra sp. ZB on Medvednica Mt. exhibited genetic distinctiveness from other samples in central Europe (Germany, Austria, Switzerland and Italy).

DNA barcoding serves as a primary tool for specimen identification and biodiversity quantification. It is increasingly utilised for the discovery and description of new and cryptic species, as well as for associating immature and mature life stages (Hebert et al. 2003a, Hebert et al. 2003b, Hebert et al. 2004). The morphological analysis in this research aligns with DNA barcoding results, confirming its effectiveness (Hebert et al. 2004, Pauls et al. 2010). The comprehensiveness of the database for target group in the researched area is crucial for linking morphologically identified individuals with DNA barcodes. Hlebec et al. (2022) published a DNA barcoding reference library for stoneflies in Croatia, facilitating precise stonefly species identification.

Our study emphasises the importance of faunistic and phylogenetic analyses for identifying populations with unique genetic traits. This underscores the necessity for increased habitat protection and conservation efforts. DNA barcoding emerges as a crucial tool in conservation biology, aiding in the development of management plans for specific animal groups and areas.

Despite its utility, DNA barcoding encounters challenges. Species delimitation may not always be straightforward, especially with variations in intraspecific genetic diversity. Some species exhibit greater intraspecific genetic diversity, as seen in Isoperla grammatica, while others show much lower intraspecific genetic diversity, such as Bracyhptera seticornis. In unexplored areas or with lesser-known organisms, the method's efficacy can suffer due to database limitations. Taxonomic experts and detailed morphological analyses are crucial in such cases. BOLD and GenBank databases are valuable resources for European stonefly fauna research, particularly in well-sequenced regions like central and northern Europe, including Medvednica Nature Park.