In total, 82 morphologically identified individuals of the W. occipitalis species complex (File S1 [Table S1]) were included in this study. Sequences for all four analysed gene fragments were available for 36 of these samples, mitochondrial information for 39 samples. All samples were collected from 1994 to 2021 in eight European countries (Fig. 2; File S1 [Table S1]). The samples were collected mostly by net sampling, occasionally with other methods such as light or malaise traps. Afterwards, the samples were temporary stored in 65 % ethanol until morphological determination and after that finally stored in 85 % ethanol at room temperature.

Morphological determination of the samples was mostly based on the morphology of the male genitalia. Some females could be determined by the co-caught males. Individuals for which no clear morphological determination to a group could be made were labelled as Wormaldia sp. Only adult individuals were included for this study.

DNA was first extracted from the whole body using the E.Z.N.A.® Tissue DNA Kit (Omega Bio-tek) following the “DNA Extraction and Purification from Tissue” protocol with minor modifications (File S1 [Table S2]). In later samples, the abdomen was removed prior to extraction because of the potential risk of damaging determination-relevant body parts and the risk of contamination from the digestive system. DNA concentration for all samples was measured with fluorescence using the DeNovix Fluorescence Assay following the DeNovix dsDNA Broad Range protocol (File S1 [Table S1]).

Four gene fragments were amplified: two fragments of the mitochondrial gene cytochrome-c-oxidase subunit I, COI5-P (658 bp) and COI3-P (551 bp), and the nuclear genes wingless (wingless, 473 bp) and nuclear rRNA and carbamoyl-phosphate synthetase 2, aspartate transcarbamylase, and dihydroorotase (CAD, 828 bp). The total volume of each PCR reaction was 15 µL. This contained 7.5 µL of Qiagen Master Mix, 1 µL (COI) to 1.5 µL (wingless and CAD) of primers, 2-4 µL of DNA and made up with water. The amplification of COI5-P was performed with the primers LCO1490 (5´-GGTCAACAAATCATAAAGATATTGG-3´) (Folmer et al. 1994) and HCO2198 (5´-TAAA CTTCAGGGTGACCAAAAAATCA-3´) (Folmer et al. 1994) or LepF1 (5’-ATTCAACCAATCAT AAAGATATTGG-3’) (Hebert et al. 2004) and LepR1 (5’-TAAACTTCTGGATGTCCAAAAAATCA-3’) (Hebert et al. 2004) using the following PCR program: 5’ 95°C, 38× (30’’ 95°C, 90” 49°C, 60” 72°C) 30’ 68°C. COI3-P was amplified using the primers Jerry (5’-CAAC ATTTATTTTGATTTTTTGG-3′) (Simon et al. 1994) and S20 (5′-GGGAAAAAGGTTAAATTTACTCC-3′) (Pauls et al. 2003) and the following PCR program: 5’ 95°C, 35× (60’’ 95°C, 90” 47°C, 90” 72°C) 30’ 68°C. The primer Wingnut1a (5’-GAA ATGCGNCARGARTGYAA-3’) (Pauls et al. 2008) and Wingnut3 (5’-ACYTCRCARCACCARTGRAA-3’) (Pauls et al. 2008) were used for the amplification of the nuclear gene wingless, 743nf-ino (5’-GGIGTIACIACIGCITGYTTYGARCC-3’) (Johanson & Malm 2010) and 1028r-ino (5’-TTRTTIGGIARYTGICCICCCAT-3’) (Johanson & Malm 2010) were used for CAD. The used PCR program of both nuclear gene fragments was 5’ 95°C, 40× (45’’ 95°C, 90” 58°C-0,2°C/Cycle, 90” 72°C) 30’ 68°C. The success of the PCR was checked by electrophoresis using 1.5 % agarose gels. PCR products were purified and prepared for sequencing using the corresponding PCR primers. The actual Sanger sequencing of the sense and antisense strands was done by Macrogen Europe (Amsterdam, The Netherlands).

The sequences were edited with GENEIOUS v9.1.8 (Dotmatics 2017), imported and manually aligned in BIOEDIT v7.2.5 (Hall 1999). Analyses were performed in R v4.2.3 (R Core Team 2023) and BEAST v2.6.7 (Bouckaert et al. 2014). The combined dataset only included samples with all four sequences available (File S1 [Table S1]). For the mitochondrial dataset, only samples with sequences available for both COI fragments were included.

For the phylogenetic analyses, first the number of haplotypes, haplotype and nucleotide diversity as well as segregating sites were calculated using the R-pages ‘pegas’ (Paradis 2010) and ‘ape’ (Paradis & Schliep 2019). Pairwise distances were calculated using TrN as the best-fitting substitution model using the R-packages ‘ape’ (Paradis & Schliep 2019) and ‘phangorn’ (Schliep 2011) and visualised in a heatmap with dendrogram. A mismatch analysis and tests of neutral mutation hypothesis with Tajima’s D for the W. occipitalis species complex, species groups and subgroups were performed using the ‘pegas’ and ‘adegenet’ packages (Jombart 2008; Paradis 2010). Furthermore, a principal component analysis (PCA) was done using the package ‘stats’ (R Core Team 2023). A TCS haplotype network was constructed using the mitochondrial dataset of the samples in R using the package ‘pegas’ (Paradis 2010).

A Maximum Likelihood tree with bootstrap values was reconstructed each for the combined, mitochondrial and nuclear dataset in R using the ‘phangorn’ package (Schliep et al. 2016). The best-suited substitution models were selected by R based on the lowest Bayesian information criterion (BIC). For the mitochondrial TrN+I, for the nuclear HKY+I and for the combined dataset HKY+G(4)+I was used as substitution models. In BEAST, a Bayesian tree of the combined dataset was constructed. The substitution model for each gene fragment was selected based on the lowest BIC calculated in JMODELTEST v2.1.9 (Darriba et al. 2012; Guindon & Gascuel 2003). The used model was HKY+G+I for COI5-P, HKY+I for COI3-P, K80 for wingless and HKY for CAD. For the molecular clock, the settings of the relaxed clock were log normal for each gene fragment, and a clock rate of 0.0177 for the mitochondrial gene fragments and 0.00177 for the nuclear genes was used (Papadopoulou et al. 2010). For the tree priors, Yule model was used for all genes. The MCMC was set to 30 million generations saving trees every 3,000 generations. The convergence and stationarity was checked with TRACER v1.7.2 (Rambaut et al. 2018). The BEAST package LOGCOMBINER v2.6.7. was used to combine the separate tree-files from the gene fragments. A consensus tree with burning percentage of 10 and posterior probability limit of 0.5 was computed using TREEANNOTATOR v2.6.7 of the BEAST package. The consensus tree was visualized in FIGTREE v.1.4.4 (Rambaut 2018). The node ages were calculated by BEAST using the same settings described for the Bayesian tree. For this purpose, two calculations per node were performed: one calculation with the samples that had the greatest and one that had the least genetic distance in the haplotype network for the respective groups. The average of these two node ages was calculated.

The genetic analyses identified two distinct monophyletic groups, W. occipitalis and W. subterranea, and four subgroups within W. occipitalis as well as three in W. subterranea (Figs 3–5). The assumed age of the split between the two major groups was estimated at about six million years (Fig. 5B). Accordingly, these two groups most likely were separated during the Miocene-Pliocene transition (Fig. 6A). During this period, climatic and geological changes occurred in the Alps, which could have promoted vicariance and speciation processes of these two groups (Barron & Keller 1982; Flower & Kennett 1994; Willett et al. 2006). For instance, heaving and depression processes that occurred in the Late Miocene in the Alps may have led to a geographical barrier in the Eastern Alps, likely resulting in range disjunction (Pfiffner 2015). A previous study also assumed diversification of the trichopteran genus Consorophylax during this time period due to the geological transformation of the Alps in interaction with climatic changes (Graf & Vitecek 2016).

During the Pliocene, these two groups were not affected by further internal differentiation. From the Pleistocene onwards, we propose several further dispersal processes causing differentiation within both species groups (Fig. 6B): On the one hand, a possible vicariance within W. occipitalis could have occurred around 650,000 years ago (Fig. 5B), initially into a south-Alpine (incl. W. occipitalis subgroups 3 and 4), and a south-west-Alpine clade (incl. W. occipitalis subgroups 1 and 2). This phenomenon is also proposed for other species in this time window (Berends et al. 2021). Each of these two groups in turn diverged into two further subgroups potentially 400‒200 kya. Thus, up to and including the Last Glacial Maximum (LGM), four geographically distinct subgroups of W. occipitalis existed, distributed across the southern/south-western Alps. On the other hand, according to the hypothesis proposed by the Bayesian tree (Fig. 5B), W. subterranea probably originally split into an Alpine and a Balkan group and afterwards spread 650,000‒470,000 years ago from the eastern to the western Alps, most likely using a northern route around the Alps. The climatic oscillation during the Pleistocene (Berends et al. 2021) apparently led to a segregation of the eastern Alpine group from the Balkan individuals and to an isolation of the western Alpine population from the eastern Alpine individuals due to glaciation.

Nevertheless, the data obtained from the molecular clock of the Bayesian analyses must be viewed critically. In contrast to the Maximum Likelihood tree, the Bayesian tree exhibits low node support. Secondly, the substitution rates established for insects were used (Papadopoulou et al. 2010). However, since Trichoptera have a lower mobility and a shorter lifespan than other insects (Curry & Baird 2015; Schröder et al. 2021) and may therefore have a different rate of molecular evolution, the applied substitution rates could be less suitable for Trichoptera.

The estimated node ages of the subgroups within both species groups call for multiple peripheral-Alpine extra-Mediterranean refugia (Fig. 6C): The W. occipitalis species groups possessed a southern central Alpine (W. occipitalis group 3), a north-western Italian (W. occipitalis group 4), a south-western Alpine (W. occipitalis group 2) and western French refugium (W. occipitalis group 1). W. subterranea apparently had three refugia in total, in the western Alps (W. subterranea group 2), eastern Alps (W. subterranea group 1) and, at least one, in the western Balkan region (W. subterranea group 3).

This biogeographic hypothesis with a number of extra-Mediterranean refugia is consistent with the phylogeographic pattern of other taxa (Malicky 2006; Pauls et al. 2006; Schmitt & Varga 2012; Schönswetter et al. 2004). Especially for semi-aquatic, stream-inhabiting taxa, such as caddisflies, extra-Mediterranean refugia, e.g. at the margins of the Alps, are more common than for terrestrial species as semi-aquatic taxa could survive glacial periods due to permanently flowing waters whose temperature never dropped below 0 °C (Malicky 1983). Such a ‘dinodal’ distribution (Malicky 1983) has been demonstrated for further (trichopteran) species (Engelhardt et al. 2011; Garcia-Raventós et al. 2021; Pauls et al. 2006). These refugia were the basis for postglacial dispersal. In our case, the western French W. occipitalis group 1 and the eastern-Alpine population W. subterranea group 1 both were spreading northwards into Central Europe postglacially (Fig. 6D; File S1 [Fig. S8]). Here, there distribution is overlapping today and even syntopic occurrences are known (Fig. 6E).

Despite not being significant, the high negative value of Taijma’s D of the W. occipitalis group 1 supports a range expansion for this group (File S1 [Table S6]); recent range expansion is also supported by the star-like haplotype structure of the mtDNA with one common central and several rare satellite haplotypes. The lack of significance of Taijma’s D might result from the limited sample size in our analysis and might have become significant if using a higher number of samples. A similar result was obtained for W. occipitalis group 3, except that Taijma’s D is strongly positive, indicating range regression (Schmitt 2020). W. occipitalis group 2 and 4 showed neither signs of range expansion nor regression. However, the sample size is not sufficient for a valid test result and should only be interpreted in conjunction with the other results.

Our genetic analyses (Figs 3–5; File S1 [Figs S6, S7]) confirm morphological differentiation and accordingly the assumption of several groups and subgroups within the species complex W. occipitalis (Neu 2015). The two morphologically distinguishable taxa W. occipitalis and W. subterranea also showed a clear genetic differentiation with a pairwise difference of almost 6 % in the combined dataset and about 10 % for the mitochondrial COI alone. Among these two taxa, the Alpine specimens do not show a stronger differentiation than the specimens from Central Europe. In general, the detected genetic distance between the two major groups of our analysis is comparable to previous studies in which an interspecific level of differentiation was assumed for members of other caddisfly families (Graf et al. 2015; Pauls et al. 2010). This high degree of differentiation might have resulted from low gene flow due to distribution in a mountain environment and generally low dispersal power of caddisflies (Curry & Baird 2015; Schröder et al. 2021). In addition, there is no detected genetic evidence of hybridisation (Fig. 3) between these two taxa since their presumed separation in the Mio-Pliocene Transition; even in the current area of sympatry, no indication for hybridisation exists, neither genetically nor morphologically. Nevertheless, potential hybridisation events cannot completely be ruled out. This strongly underlines that reproductive isolation may also have fostered the high degree of genetic differentiation between these two groups that definitively should be considered as good species: W. occipitalis and W. subterranea.

In contrast to morphology, the genetic analyses revealed four instead of two subgroups within W. occipitalis and additionally a differentiation into three subgroups within W. subterranea (Figs 3–5). The morphologically determined subspecies W. o. occipitalis and W. o. meridionalis were also identified as groups. However, W. o. meridionalis did not represent a monophyletic group but was observed in three genetic subgroups including two monophyletic clades. The morphological group determined as hybrid W. occipitalis x subterranea was genetically not detected to be of hybrid origin. Instead, the specimens were clearly classified within W. occipitalis, forming a monophyletic clade there. The grouping of individuals was mostly consistent between the genetic and morphological determination (Figs 4, 5). Only the specimens TRI0063 and 64 (France, Chamrousse) morphologically grouped within W. o. meridionalis but were assigned to the W. o. occipitalis cluster by genetic analyses. All specimens morphologically not determined to the species level (Wormaldia sp.) clearly clustered within one of the genetic groups; hence no further lineages were recovered. The mitochondrial genetic distance among the subgroups was 1.4‒2.8 % within W. occipitalis and about 3.0‒3.5 % for the three subgroups of W. subterranea. Similar intraspecific genetic distances of genetic clusters were previously found within other trichopteran families (Graf et al. 2015; Previšić et al. 2014; Previšić et al. 2009). Higher intraspecific genetic distances have also been detected (Pauls et al. 2006), however, here species status needs to be re-evaluated (Oláh et al. 2015).

Geographic proximity partly correlated with genetic differentiation and phylogenetic relatedness of the species’ subgroups: On the one hand, the subgroups of W. occipitalis show an exact order in the haplotype network (Fig. 4) according to their geographical distribution from the western to the southern Alps. Hence, they fulfil all criteria of allopatric subspecies. On the other hand, the German specimens within W. subterranea group 1 are geographically more distant from the Austrian specimens of the same phylogenetic group than from the genetically more differentiated specimens of W. subterranea group 2. Nevertheless, their subspecific status, as well as for W. subterranea group 3, is highly likely.

In short, the taxa W. occipitalis and W. subterranea are distributed in the same area and habitat in western Germany (Fig. 6E). This syntopic distribution is an indication that these two taxa are distinct species, as previous morphological studies already assumed (Neu 2015; Oláh et al. 2019). The high degree of genetic differentiation, lack of genetic evidence of hybridisation and the morphological distinctness support this hypothesis. The subgroups of W. occipitalis and W. subterranea are genetically distinct and have allopatric distributions. However, difficulties in the unambiguous morphological assignment of individuals to a subgroup as well as poor bootstrap values or Bayesian posterior probabilities in phylogenetic analyses and relatively low genetic distances tend to argue for a differentiation at subspecific level.