Adults of I. popijaci sp. nov. were collected in June 2019 at the entrance to the Ševerova Cave (karstic source of the intermittent Krasulja rivulet in Krbava field). A subsequent collecting trip upstream of the Krasulja rivulet (in June 2021), near the karstic source adjacent to the village of Mirići, resulted in finding more specimens of I. popijaci sp. nov.

A total of 42 specimens (34 adults and 8 larvae) belonging to Isoperla popijaci sp. nov., were collected. Adult specimens were collected using sweep nets, while larval specimens were collected by handpicking. The aedeagus was everted in the field and specimens were fixed and stored in 96% ethanol for morphological and molecular analysis. Morphological characteristics of male terminalia were examined after potassium hydroxide (KOH) treatment.

The holotype and part of the paratypes series are deposited in the Croatian Natural History Museum, Zagreb, Croatia (CNHM), Collection of Plecoptera Sivec & Hlebec, while other paratypes are kept in the Slovenian Museum of Natural History, Ljubljana, Slovenia (PMSL).

Photographs, diagnostic characterisation and comparative morphological examination of specimens were made using a ZEISS SteREO DiscoveryV.20 stereomicroscope. Pencil drawings were produced with a camera lucida and then digitally edited and inked. Figures 3A, B, 4A–D (SEM images) were made using a JEOL JSM-7000F scanning electron microscope. The penis (one of paratype specimen) for the SEM study was critical-point dried (Figure 4A–D).

Nomenclature is in accordance with the International Code of the Zoological Nomenclature (ICZN 1999). The species is proposed by following the rules of the Code. Abbreviations for the type specimens are HT–holotype, PT–paratype and PTs–paratypes.

Comparative study on the morphology of penial structures was conducted using ten species belonging to the genus Isoperla, collected in Croatia: I. bosnica Aubert, 1964; I. inermis Kaćanski et Zwick, 1970; I. rivulorum (Pictet, 1841); I. lugens (Klapálek, 1923); I. illyrica Tabacaru, 1971; I. tripartita Illies, 1954; I. grammatica (Poda, 1761); I. difformis (Klapálek, 1909); I. oxylepis (Despax, 1936) and I. albanica Aubert, 1964. Morphological taxonomic classifications follow the traditional system (Poda 1761; Pictet 1841; Klapálek 1909, 1923; Despax 1936; Illies 1952, 1954, 1966; Aubert 1964; Tabacaru 1971; Kaćanski and Zwick 1970, Murányi 2011; Murányi et al. 2016).

One male, one female and one larva of Isoperla popijaci sp. nov. were used in molecular analyses and mutually associated. DNA was extracted from the single leg of specimens using QIAamp DNA Micro Kit (Qiagen, Germany) according to the manufacturer’s specifications and eluted in 50 µl of elution buffer. The 5’ fragment of the mitochondrial cytochrome c oxidase subunit I gene (COI) was amplified using standard PCR-protocols and four sets of primers: LCO-1490/HCO-2198 (Folmer et al. 1994) or C_LepFolF/C_LepFolR (as was used in Hebert et al. 2004) or a combination of MLepF1/LepR1 and MLepR1/LepF1 (yielding two shorter, overlapping fragments as was used in Hajibabaei et al. 2006) in 20 µl reactions. Polymerase chain reactions (PCRs) for all primer sets were carried out using: 1 x DreamTaq reaction buffer with 2 mM MgCl₂ (Thermo Fisher Scientific Inc., US), 0.2 mM dNTPs, 0.4 µM of each primer, 0.025 U/µl of DreamTaq polymerase (Thermo Fisher Scientific Inc., US) and 1 µl of eluted DNA. For the first mentioned primers set (LCO-1490/HCO-2198) the following PCR cycling conditions were applied: initial denaturation at 95°C for 2 min, followed by 35 cycles of denaturation at 95°C for 30 s, annealing at 50°C for 30 s, extension at 72°C for 1 min, followed by a final extension step at 72°C for 10 min. PCR products were purified using Exonuclease I (0.05 U/µl), FastAP Thermosensitive Alkaline Phosphatase (0.025 U/µl) enzymatic system (Thermo Fisher Scientific Inc., US). The reaction was carried using the protocol: 1 h at 37°C followed by 20 min at 80°C. Sequencing was performed by Macrogen Inc. (Amsterdam, The Netherlands) using the amplification primers. Sequences obtained in the study were deposited in the BOLD database (Ratnasingham and Hebert 2007) and GenBank (under the accession numbers MW907977–MW907980, MW907982–MW907988 and MW907990–MW907993).

In total, 15 obtained Isoperla sequences were checked, edited, assembled from both directions and inspected manually for base-pair ambiguities, as well as stop codons, indels or double peaks in chromatograms (as indicators for the possible erroneous amplification of nuclear mitochondrial pseudogene) in Geneious R6 (https://www.geneious.com). All available Isoperla sequences were retrieved from the GenBank and BOLD databases (accessed 10/01/2021) and aligned with sequences from this study using MAFFT v.7 (Katoh and Standley 2013). Any length variants were excluded from the final alignments. Sequences were collapsed into 456 unique COI haplotypes using the online tool FaBox v.1.5 (Villesen 2007) and, from all species, the most diverse haplotypes from I. tripartita and I. rivulorum species group, as well as species I. lugens, were retained for further analysis. The final dataset for phylogenetic analysis and species delineation comprised 27 sequences, including 10 haplotypes observed in this study (see Table 1). Isoperla obscura (INTAP055-17) and Taeniopteryx burksi (08INHSP-002) were selected as outgroups according to the North American Plecoptera phylogeny published by South et al. (2020). Amongst morphologically-defined species, evolutionary divergence was estimated using the pairwise comparison of the uncorrected genetic distances (p-distances) in MEGA-X (Kumar et al. 2018). For p-distances, a colour heat map was drawn using the Python data visualisation library Seaborn (version 0.11.1, Waskom 2021). Phylogenetic relationships were estimated by three different optimality criteria: Neighbour Joining (NJ), Maximum Likelihood (ML) and Bayesian Inference (BI). NJ and ML were performed in MEGA-X (Kumar et al. 2018), while BI in MrBayes 3.2.7. (Ronquist et al. 2012). For ML and BI, the optimal model of nucleotide evolution (Hasegawa-Kishino-Yano model with gamma distributed rate variation amongst sites and a significant proportion of invariable sites: HKY+I+G) was selected under the Bayesian Information Criterion (BIC) using jModelTest 2.1.5 (Darriba et al. 2012). Nodes in the phylogenetic trees with bootstrap values P ≥ 70 in NJ and ML and posterior probabilities values pp ≥ 0.90 in BI were considered well supported. NJ was made using the Kimura-2-parameter (K2P) model of nucleotide substitution with the pairwise deletion option. Bootstrap support was inferred using the fast bootstrap algorithm, based on 5000 replicates. Nearest-Neighbour-Interchange (NNI), a heuristic method using the fast bootstrap algorithm, was used in ML with 1000 replicates.

For BI, the dataset was partitioned by codon positions. Two separate runs with four Metropolis-coupled Monte Carlo Markov chains (MMCM) were performed for 10 million generations while trees were sampled every 1000 generations with the first 25% of sampled trees discarded as burn-in. The remaining trees were used to create a 50% majority rule consensus tree. TRACER v.1.7.1 (Rambaut et al. 2018) was used to check the convergence between the two runs. The phylogenetic trees were visualised using FigTree v.1.4.3. (Rambaut 2009) and iTOL v.5 (Letunic and Bork 2021). Several methods of species delimitation were applied: the Automatic Barcode Gap Discovery (ABGD) method (Puillandre et al. 2012), the Bayesian implementation of the Poisson Tree Processes (bPTP) method (Zhang et al. 2013) and the multi-rate Poisson Tree Process (mPTP) method (Kapli et al. 2017). The ABGD was performed at the web server by using the K2P model. All values were set to default, except the value of relative gap width, which was set to 1, while the default gap width of 1.5 resulted in a single group. The bPTP method was performed on the web server at http://species.h-its.org, while the mPTP method was run on the web server at http://mptp.h-its.org/. Both methods were applied using default parameters, outgroups have been removed from the analysis and the same ML input tree was used.

The new species should probably be regarded as Critically Endangered (CR) or Vulnerable (VU) by the IUCN Criteria. Up to now, it is known only from the areas nearby two karstic sources.

The alignment of COI gene sequences was 658 bp in length and comprised of 202 variable sites, of which 139 were parsimony informative. Three implemented criteria of phylogenetic reconstruction (NJ, ML and BI) resulted in congruent topologies with highly similar support values (Figure 6), characterised by the presence of two deeply divergent lineages, I. popijaci sp. nov. and “Isoperla PL”, which did not cluster with any of the currently defined taxa.

Figure 6. 

Maximum Likelihood cladogram, based on the analysis of the COI haplotypes of Isoperla species. Numbers at the nodes indicate Neighbour-Joining (NJ), Maximum Likelihood (ML) bootstrap support values (BS) and Bayesian posterior probabilities (BPP), respectively. The results of species delimitations are represented with the vertical bars, from left to right, indicate the OTUs inferred by bPTP, ABGD and mPTP. “Isoperla PL” indicates additional separate lineage obtained in this study. Terminal codes present BOLD/GenBank Process ID, as in Table 1.

Mitochondrial COI sequences, obtained from I. popijaci sp. nov. (adults and larva), were identical (a single unique haplotype). The monophyly of the newly-described species is highly supported (Figure 6). This species represents the first branch-off within the clade comprised of monophyletic I. lugens and I. rivulorum subclades, as well as another tentative new taxon obtained in this study (clade designated as “Isoperla PL” with representatives CROPL214–21 and CROPL230–21). The designation “PL” denotes the abbreviation Plitvice Lakes, nearby where a specimen was found. Five sequences of “Isoperla PL” represent 2 haplotypes (CROPL214–21 and CROPL230–21) with low intraspecific uncorrected p-distance (0.0096).

Intraspecific uncorrected p-distances are as follows for the following species: 0.32–1.59% in I. rivulorum, 0.16–0.48% in I. lugens, 0.01–7.82% in I. tripartita, 0.32% in I. vjosae and 0.16% in I. illyrica. Interspecific uncorrected p-distances for I. popijaci sp. nov. ranged from 6.69–12.59%; specifically, 6.69–7.17% to I. rivulorum, 8.15–8.45% to I. lugens, 9.99–10.22% to I. vjosae, 10.4–12.6% to I. tripartita, 10.38–12.61% to I. illyrica, 10.69% to I. pesici and 8.12% to the “Isoperla PL” (Figure 7). Overall, observed intraspecific genetic distances within the genus ranged from 0.01–7.82%.

Figure 7. 

Colour heat map showing inter- and intraspecific uncorrected p-distances of the mitochondrial cytochrome oxidase subunit I (COI) barcode region. Isoperla popijaci sp. nov. and Isoperla PL appear as highly divergent. Intraspecific p-distances are outlined by the black line.

Within the I. rivulorum clade, Croatian sample CROPL066–21 appeared as a separate lineage, subdivided from Alpine specimens (Figure 6). The uncorrected p-distances for sample CROPL066–21 are in the range 1.28–1.59% to other I. rivulorum samples.

A well-supported clade comprised two newly-discovered lineages (Isoperla popijaci sp. nov. and “Isoperla PL”), together with I. lugens and I. rivulorum, and was recovered in all three tree-building algorithms.

According to the results of the first molecular characterisation of I. illyrica obtained in this study, specimens clustered in a within the monophyletic clade with intraspecific uncorrected p-distances of 0.16%. Interspecific p-distances between I. illyrica and I. tripartita ranged from 0.96–5.91%.

All species delimitation analyses (bPTP, ABGD and mPTP) for mtDNA (COI) have delineated two well-separated lineages Isoperla popijaci sp. nov. and “Isoperla PL” as tentative species. Applied methods resulted in various numbers of delineated groups. In the ABGD analysis, initial partitioning identified eight, while recursive partitioning showed the existence of nine putative species for the majority of prior intraspecific divergence values (P). The mPTP method delimited seven operational taxonomic units (OTUs) and, according to these results, is the most conservative approach, while the bPTP recognised 9 OTUs.

Contrary to ABGD and bPTP, the mPTP analysis shows I. illyrica, I. tripartita and I. pesici, morphologically assigned to I. (tripartita) species group, as a single OTU. These species are completely separated into three OTUs in the bPTP analysis. The separation of sample CROPL122-21 as a distinct species (I. tripartita) was supported by all three species delimitation methods.

The lowest interspecific p-distance between I. popijaci sp. nov. and I. rivulorum was found to be 6.69%, indicating distinct species. This exceeds intraspecific divergences (ISD ≥ 2%) commonly used as one of the criteria for a delimitation of closely-related species in aquatic insects: Ephemeroptera, Plecoptera and Trichoptera (Ball et al. 2009; Zhou et al. 2009). Values above 2% have already been reported amongst Plecoptera (Zhou et al. 2010, Gill et al. 2015), which was probably caused by poor mobility of some Plecoptera species (Boumans and Baumann 2012) and, consequently, geographical isolation among populations.

The finding of the second well-separated lineage (“Isoperla PL”), most closely related to species I. rivulorum (interspecific p-distance from 6.54–7.19%) implies existence of another new species of the genus Isoperla (unpublished data). Taxa obtained in this study (Isoperla popijaci sp. nov. and “Isoperla PL”) are separated by a large interspecific p-distance of 8.12%. Future research will seek to determine whether this value has repercussions to the geographical isolation and specificity of the (micro-) habitats in which the taxa were found.

Based on the occurrence of I. lugens (alpine species) and I. rivulorum (alpine, central European species) in the Dinaric karst and their appearance as the most recently diverged lineages within I. popijaci + “Isoperla PL” + I. lugens + I. rivulorum clade (Figure 6), it can be assumed that the Dinaric karst might represent the area of origin of those alpine species as well as the diversification centre from where they spread northwards. However, to test this hypothesis, data across the whole distributional range and use of other molecular markers (mitochondrial and nuclear as well) are necessary.

To establish a final phylogenetic relationship in the monophyletic I. tripartita species group, it is necessary to collect specimens from its entire range and use a multi-gene molecular approach as well.

Previous research showed the wide range of variability in intraspecific divergence within the order Plecoptera (Zhou et al. 2009; Gill et al. 2015; Stark et al. 2015) and uncorrected intraspecific p-distances from our study (0.01–7.82%) are consistent with the previously reported values.

Based on the morphological characteristics, the new species can be assigned to the Isoperla tripartita species group. The I. tripartita species group is characterised by the divided medial penial armature (into upper and lower coloured portions, divided or subdivided) and lateral penial armatures (Illies 1954; Murányi 2011; Murányi et al. 2016). Popijač’s Yellow Sally is characterised by divided medial penial armature, with the distal part bearing short spines, but with indistinct lateral penial armature. The genetic distinction, in combination with morphological features, is significantly different from all other species and promotes I. popijaci sp. nov. as a new species.

Phylogenetic reconstructions support the monophyly of the I. tripartita species group, which is, together with I. grammatica, notable by the high morphological variability of certain species (Zwick 1978; Murányi 2011). In Croatia, significant morphological variability has been also observed in I. inermis from different localities (personal observation), of which some are very similar to I. difformis (Central European species) in the penial armature. Therefore, future studies should investigate relationships between and within Isoperla populations from the Balkan Peninsula (e.g. Cetina River, National Park Plitvice Lakes, Kupa River and nearby springs in Slovenia) by applying a multi-gene approach.

Other species are somewhat less variable and occupy smaller distributional areas (as recently described species from Europe and Asia). Those endemics are of special interest to our study because it is assumed that more endemics species are likely to be discovered, especially in poorly-explored areas with high biodiversity like the Balkans. More new species are expected to be found in Croatia, as the majority of the country’s territory has not been studied yet regarding Plecoptera.

Anthropogenic activities have already resulted in the reduction of population size (especially larger species from the genera Perla, Dinocras and Perlodes) (personal observation). All the above-mentioned calls for more detailed studies of species distributional patterns, as well as of genetic diversity of populations. Emphasis should also be put on the isolated habitats (karst areas) as they can have the highest conservation value as refugium and the maintenance of genetic diversity.

Until now, Popijač’s Yellow Sally is known to inhabit the parts of the rivulet close to two karstic sources, of which one is a cave entrance. Although there are no true troglobionts within the order Plecoptera, several species have been found to inhabit stream sources around the openings of caves (for example I. inermis) and there are no records of these species from the downstream part of the same stream. Another example is Brachyptera tristis (Klapálek, 1901), a species that spends its entire life cycle underground (the stream of Krupa River) (personal observations). It is, hence, important to pay special attention to the research of caves, pits, underground and temporary rivers and streams that abound in the Dinaric karst geology. These habitats host some of the most complex and diverse faunas (Culver and Sket 2000) as a consequence of composite geological history and the intensive process of karstification (Sket 1999). The Balkan Peninsula is known for its high biodiversity (Sket et al. 2004), especially of aquatic species (Kryštufek et al. 2007, Previšić et al. 2009, 2014; Murányi 2011; Vitecek et al. 2015; Kučinić et al. 2017). It can be expected that future research will contribute to the discovery of biodiversity patterns as well as new species, especially microendemic species (Graf et al. 2009, 2012; Kučinić et al. 2013; Vitecek et al. 2017). Karst habitats, such as Ševerova Cave, represent some of the most dynamic freshwater habitats, especially in terms of biological-geological interactions (Ridl et al. 2018). With the alternation of wet and dry phases and temporal dynamics of water flow, temporary rivers have a great influence on local ecological interactions, both in aquatic and terrestrial habitats (Larned et al. 2010). It is a significant assumption that climate change will increase the duration and frequency of dry phases, so it is expected that this will lead to the disappearance of taxa whose entire life cycle (or at least part of it) is related to aquatic environments (Larned et al. 2010).