Specimens were collected with light traps or sweep nets and preserved in either 70% or 95% ethanol. Prior to DNA extraction, wings (in adults), heads and terminal abdominal segments were removed as specimen vouchers. DNA was extracted with hot sodium hydroxide as described in Truett et al. (2000). Extracted DNA was visualized on an agarose gel to estimate the mean size of the fragments. Concentration of double stranded DNA was flourometrically estimated using a Qubit assay kit (Thermo Fisher Scientific). The DNA pellet was washed 2x with 70% EtOH. Following the wash, the EtOH was discarded and the pellet dried. The pellet was sent dried to the Center for Anchored Phylogenomics at Florida State University (http://www.anchoredphylogeny.com) for quality check and library preparation, sequencing and raw data handling following Lemmon et al. (2012) and Prum et al. (2015).

The probes were designed in collaboration with the Center for Anchored Phylogenomics following Breinholt et al. (2018) and Deng et al. (2022). In short, the transcripts identified as homologous with the existing Lepidoptera probe set from the 1Kite Project (Misof et al. 2014, Breinhold et al. 2018) were aligned in MAFFT v7.023b (Katoh & Standley 2013) by target locus, then trimmed to conserved regions (Lemmon et al. 2012), and finally manually inspected in Geneious R9, (Biomatters Ltd., Kearse et al. 2012). A total of 960 target loci (averaging 232 bp in length) remained after masking/removing regions identified to be repetitive using k-mer distribution profiling (see Hamilton et al. 2016 for details). For each of these loci, probes were tiled uniformly across each of the sixteen reference sequences with 4.2x coverage, to produce 57094 probes. Probes were produced by Agilent as a SureSelect XP kit. Methods for probe tiling are described in Deng et al. (2022).

Extracted DNA was sonicated using the Covaris Ultrasonicator to a size distribution of 150–400 bp. Libraries were then prepared following Lemmon et al. (2012) and indexed with single 8-bp indexes chosen for being different at two or more sites. Libraries were then combined in pools of 16 samples and enriched using the Agilent® SureSelect XP probe kit. The enriched libraries were sequenced on two lanes of an Illiumina HiSeq 2500 sequencer with a paired end 150 bp sequencing protocol and onboard cluster generation. The raw data totaled 88 gbp.

Raw read pairs passing the CASAVA high-chastity filter were merged following Rokyta et al. (2012) and library adapters were trimmed using TrimGalore! V. 0.4.0 (https://github.com/FelixKrueger/TrimGalore). Reads were then assembled to the reference R. fasciata (probe region sequences) using the assembly procedure described by Prum et al. (2015) with the same parameters as Hamilton et al. (2016). Consensus sequences were derived from assembly clusters containing ten or more reads and ambiguities were called when site patterns could not be explained by a 1% sequencing error. Orthology among consensus sequences was determined using a neighbor-joining approach based on a pairwise distance matrix (see Hamilton et al. 2016 for details). Orthologous sequences were aligned using MAFFT v7.023b (Katoh & Standley 2013) in Geneious R9 (Biomatters Ltd., Kearse et al. 2012).

High-throughput sequencing data have been shown to be susceptible to contamination (e.g.Lusk 2014). Also, poor alignment of individual taxa or entire regions within aligned loci impact the phylogenetic inference derived from those alignments. We therefore performed multiple filtering steps to obtain our final data set for analysis.

We used Gblocks (Castresana 2000) to remove poorly aligned regions of regions containing gaps as we expected invariant length and gap free alignment across the exons we sequenced. We used default settings to obtain a stringent alignment without indels.

We then used a relative branch length approach to identify potential contaminants (Fig. 2). The rationale for this is that if a taxon is constrained in a phylogeny to a certain position, its relative branch length in each locus should vary following a normal distribution. In the case that a sequence represents a contaminant and is forced into a phylogenetic position it would not take based on the data, that branch should be unusually long. Thus, the branch will be an outlier to the normal distribution and could be interpreted as the result of contamination. We used the outlier definition according to Tukey (1997): it uses the first and third quartile (Q1, Q3) of the branch length distribution and the interquartile range IQR = Q3 – Q1. An outlier is defined as a value that is larger than Q3 + IQR * k or smaller than Q1 – IQR * k, with k = 1.5. Because this criterion is very conservative and can lead to false positives, we used k = 3. The implementation also allows us to flag only extreme outliers with k = 9, and the option to use a flexible σ multiplier instead, with outliers being outside of [Q1 – k * σ, Q3 + k * σ].

Prior to using this approach on our study data, we examined its utility and reliability on simulated DNA sequence data. First, we used Seq-Gen v.1.3.3 (Rambaut & Cassidy 1997) to simulate DNA sequences for 100 loci (500 bp, HKY model, nucleotide frequencies A: 0.294, T: 0.246, G: 0.241, C: 0.221) based on the characteristics of our own concatenated alignment and the resulting guide tree. To simulate contamination, we manipulated the simulated data set in three ways: 1) by randomly replacing 10 complete loci with others from within the simulated data set to imitate cross contamination within a study; 2) by randomly inserting 10 120-bp-fragments from within the simulated data to other species sequences to imitate chimera production in the original sequence assembly; and 3) by randomly inserting 10 120-bp-long fragments of random data to imitate contamination from an unknown source in the original sequence assembly. Randomization of loci, taxa and insert position were based on Random Thing Picker (http://andrew.hedges.name/experiments/random/pickone.html, accessed on May 12, 2016). We then performed phylogenetic analysis on all 4 simulated data sets (each form of contamination individually and all contaminations combined) under the topological constraint and estimated relative terminal branch lengths for each taxon individually across all loci. We used the approach described above to identify if the integrated contaminants resulted in outlier branch lengths. Under the topological constraint of our guide tree, our approach identified 25 of the 30 contaminations (all random inserts, all cross contaminations, five of the 10 chimeras) from our simulated data set. We also flagged two false positive outliers. These results show that the approach is useful in detecting contaminants, including small random inserts and potentially even chimeras with congeners. However, half of the chimeras were not flagged and it is possible that some false positives could be flagged and uncritically removed from a data set. We view the approach as a tool for flagging potential contaminants that allows us to further scrutinize the flagged sequences by hand, since it is time-consuming and error-prone to scrutinize all sequences of all loci manually.

We applied our contamination detection approach to our original 61 taxon data set. For the input constraint tree, we used a maximum likelihood tree generated in RAxML v. 8 (Stamatakis 2014) using the concatenated data set under a GTRCAT model with 100 fast bootstrap replicates. We then calculated the distribution of relative maximum likelihood branch lengths of each taxon for all 451 loci and flagged all sequences that represented statistical outliers in relative branch length for each taxon. The procedure was straightforward for the 66 loci that had sequence data for all 61 taxa. For 385 loci we had some taxa with missing data (from 1 to 23 taxa missing). For each locus with missing data, we generated a new constraint tree by dropping the terminals for which we had no data and used the locus-specific reduced constraint trees to identify potential contaminants using the method described above.

Our approach flagged 33 sequences from the 66 loci that had sequence data for all 61 taxa, and 230 sequences from the 385 loci with missing data, resulting in 263 total flagged sequences. We removed these sequences unless the outlier branch length was the result of the sister taxon to the flagged taxon having been removed from the single locus constraint tree due to missing data. In this case, the terminal branch length of the flagged specimen is necessarily longer than usual because the closest taxon is missing from the single locus constraint tree. Using these guidelines we removed 126 potential contaminant sequences.

Additionally, we removed 116 potential cross-contaminated sequences that were identical among non-sister taxa (based on the constraint tree). In total 242 individual sequences were removed (1.03%).

The most important predictor of phylogenetic informativeness of a particular locus is simply sequence length (Kjer et al. 2007). We thus concatenated remaining loci into “super loci”, i.e., we concatenated all loci that mapped to the same gene, but different exons in Bombyx mori. Following this, we removed all (super) loci that were <100 bp in length as well as all (super) loci and taxa with >35% missing data.

This filtering resulted in a final data set comprising 60 taxa, 159 (super) loci ranging from 302-1922 bp in length (Ø 570.5 bp), and a total concatenated alignment length of 90,717 bp. The molecular data sets are available at the Senckenberg (Meta) Data Portal under https://dataportal.senckenberg.de/dataset/cf47cdbb-7384-47e4-a89d-cb1abcf65963.

We followed three general approaches to infer the phylogenetic relationships in our data set. First, we generated a maximum likelihood tree on the concatenated alignment (90,717 bp). We used the PartitionFinder 2.0 relaxed clustering algorithm with default settings to select an optimal partitioning scheme (Lanfear et al. 2017). Then we selected nucleotide substitution models for each subset of the optimal partitioning scheme using ModelFinder as implemented in IQtree v. 1.6 (Kalyaanamoorthy et al. 2017, Nguyen et al. 2015). We then used to IQtree v. 1.6 to estimate 50 maximum likelihood trees, 25 replicates with a parsimony tree and 25 with a random starting tree, and chose the topology with the best likelihood score. We evaluated node support using 1000 ultrafast bootstrap replicates in IQtree v1.6 (Hoang et al. 2018).

Second, we generated a species tree in ASTRAL v. 5.6.1 (Zhang et al. 2018), which accounts for incomplete lineage sorting by considering conflict among gene trees. For each (super) locus, we selected the best fit model in IQtree using ModelFinder, and then generated 10 maximum likelihood trees and chose the topology with the best maximum likelihood score. We then used the best topology for each (super) locus as input into ASTRAL to generate the species tree (Supplementary Material B).

Third, we identified the best partitioning scheme and nucleotide substitution models via ModelFinder. For each of the 159 (super)loci, we used independent, variable substitution rates assuming an average substitution rate of 1 across all (super)loci and independent character state frequencies, transition/transversion rate ratios, gamma parameter shapes and proportions of invariable characters. Based on this model we generated a Bayesian tree with MrBayes v. 3.1 (Ronquist & Huelsenbeck 2003). We ran Bayesian inference for 10 x 10⁶ generations in 2 independent runs with 8 chains each, sampling each 5000^th generation. We used average standard deviations of split frequencies as reported by MrBayes and visualizations of MrBayes parameter files in combination with ESS estimates as available in Tracer v. 1.6 (Rambaut et al. 2018) to determine when analyses had converged and reached stationarity. After roughly 4.5 × 10⁶ generations, average deviation of split frequencies remained under 0.01 while visualizing MrBayes parameters indicated stationary distribution with ESS estimates generally above 200 for the majority of parameters when using a burn-in of 50%. Accordingly, we discarded the first 1000 trees as burn-in and used the remaining trees for ancestral character state estimations (see below) as well as to construct a strict consensus tree.

In total, we differentiated the eleven known gill types in Rhyacophilidae that were included in our phylogeny. Species where larval gill types are not known were coded as “unknown” and treated as missing data. We estimated the ancestral character state of gill type at each supported node along the backbone of our phylogeny as well as all nodes within the R. vulgaris group. Ancestral character estimation (ACE) was performed in BayesTraits V3.0 (available from http://www.evolution.rdg.ac.uk), using the BayesMultistate method (Pagel et al. 2004). The multistate method estimates trait values at ancestral nodes for traits that adopt a finite number of discrete states using Markov chain Monte Carlo (MCMC) or Maximum Likelihood (ML) methods. This method fits a continuous-time Markov model to the discrete character data. In this model, trait values can change along the branches at infinitesimally small time intervals, and transition rates from state i to j are defined by a transition rate q_ij.

We performed MCMC ancestral character estimation on 16 nodes with a uniform prior on a sample of trees from the MrBayes output at two levels: first with a binary coding of gills/no gills, and second for all eleven character states (Figs 1, 3).

Abdominal tracheal gills in caddisflies have been associated with osmoregulatory and breathing functions (Badcock et al. 1987). We thus assessed the interaction between gill types and the river zonation preference of larval Rhyacophila for the European species only. Taking advantage of the longitudinal river zonation as a general proxy for oxygen, current, and conductivity conditions (Moog et al. 2002), Ofenböck et al. (2010) developed the longitudinal river zonation index (LZI) to categorize European benthic invertebrate species according to their longitudinal occurrence within European river systems. The ecological analysis was necessarily restricted to European species of Rhyacophila because detailed species-level information on gill types and ecological traits are not available or known in other regions. We used the data from freshwaterecology.info (Graf et al. 2015, Schmidt-Kloiber & Hering 2015) for this analysis. Of the 116 species listed, the gill types and river zonation are known for 50 species (Supplementary Material C). We applied analysis of variance (ANOVA) to assess if river zonation occurrence differed between the different gill types, the null hypothesis being that all gill groups have the same river zonation index distribution.

Our final data set comprised 60 taxa (Supplementary Material A). Of these, 19 species belonged to the vulgaris-species group sensuSchmid (1970). Thirty-nine additional species of Rhyacophila were sampled from all evolutionary “branches” proposed by Schmid (1970) covering 11, 5, and 10 additional species groups from North America, Europe, and Asia, respectively. In addition, our data set included one specimen of Himalopsyche kuldschensis (Rhyacophilidae) and one Glossosomatidae (Agapetus hessi) as outgroups. Within our taxon sampling the gill configuration of last instar larvae is known for 43 species (~72%). In our analyses we included all known gill configurations.

Briefly, all phylogenetic analyses recovered almost congruent tree topologies: only minor changes in branch length were returned and almost all relationships received high support values (Fig. 3). We thus describe the results on the basis of the Bayesian Inference, which was also used for ancestral character state reconstruction. Phylogenetic groupings, i.e., ‘branches’, ‘veins’ and ‘species groups’ were introduced by Schmid (1970) as hierarchical categories with branches forming large clades, veins forming groups within branches, and species groups forming smaller groups within veins. Unless stated otherwise, we use this terminology, and compare our results with the branches and species groups hypothesized by Schmid (1970). Regarding the deep splits within Rhyacophilidae, R. rickeri is sister to all other Rhyacophila and the nested Himalopsyche kuldschensis. The Rhyacophila stigmatica group follows as sister to all remaining Rhyacophila. At this point the phylogeny splits Rhyacophila into two main clades. Clade A (pp = 1.0) comprises species from the divarica, naviculata, philopotamoides and vulgaris branches, and Clade B (pp = 1.0) comprises species from the naviculata, philopotamoides and vulgaris branches, rendering these branches as defined by Schmid (1970) polyphyletic. The deeper nodes in clades A and B resolve the relationships within these clades but exhibit very short internodes. At the species level, six of seven species groups sampled with more than one species were inferred as being monophyletic (except the naviculata group).

Of the vulgaris branch species in Clade A, R. oreata, R. dongkyapa and R. petersorum form a monophylum. The remaining vulgaris branch species form a clade (Clade C; pp = 1.0) within Clade B. In Clade C, R. fuscula of the fuscula group is sister to the supported vulgaris species group clade. Within the vulgaris group (Clade D-G; all pp = 1.0) all relationships are resolved. However, some internodes are very short.

Both structurally simple and complex gill types are distributed across the phylogeny, with only a few gill types representing apomorphies for specific branches, veins or species groups: H. kuldschensis, R. alberta, R. bonaparti and R. coloradensis have unique gill types in the phylogeny, while gill types observed in R. alberta, R. coloradensis and H. kuldschensis are also known from other species that we did not include in our study (Giersch 2002, Hjalmarsson et al. 2018). Single filament gills (ProsrhyacophilasensuDöhler (1950)) occur in different clades within Clade B. Gill tufts occur in several taxa of Clade B and can be differentiated depending on the configuration in which they occur on each abdominal segment: in R. brunnea (brunnea group sensuSmith 1995) there are three independent tufts, in R. sequoia (grandis group sensuSmith 1995) there are two independent tufts and the vulgaris group exhibits a single tuft (Rhyacophilasensu stricto as defined by Döhler 1950). Four digitate gills (Pararhyacophila as defined by Döhler (1950) and comb-like gills (Hyperrhyacophila as defined by Döhler (1950)) are limited to the vulgaris group.

In the vulgaris group (Fig. 3, clade D), R. loxias (single filament gill) is sister to all other taxa. Rhyacophila ferox subsequently follows as sister to the remaining vulgaris group species. Rhyacophila ferox has four digitate gills on each abdominal segment. In our analysis it is placed within the vulgaris group, though its placement was previously considered to be isolated or distantly related to R. fasciata on the basis of preanal appendages, anal sclerites and the aedeagus (Graf 2006; Graf et al. 2009).

The vulgaris group then splits into three clades that correspond to different gill configurations. The Pararhyacophila clade (Clade E) comprises R. rectispina, R. intermedia, R. lusitanica and R. munda, whose larvae exhibit four digitate gills on each abdominal segment. The Pararhyacophila clade is sister to the Hyperrhyacophila + sensu stricto clades. The Hyperrhyacophila clade (Clade F) comprises all sampled species with larvae exhibiting comb-like gills (R. torrentium, R. evoluta, R. occidentalis) as well as the sensu stricto species R. moscaryi at its base. All other species with tufted gills in the data set fall into the sensu stricto clade (Clade G). Based on the polarization of the tree with the outgroup, it appears that the clades within the vulgaris group follow an evolutionary trajectory where deeper nodes exhibit ancestors with simple gill configurations and shallow nodes bear ancestors with complex gills.

The IQTree ML analysis resulted in a topology that showed the same topology as the Bayesian Inference. We highlight differences in support values on those nodes with bootstrap supports (bs) lower than 70 (Fig. 3). The ASTRAL species tree of the 159 loci (each > 300 bp in length; Supplementary Material B), shows almost identical relationships throughout. The only conflict of supported nodes (bs > 70) with the Bayesian and ML inference relates to the placement of H. kuldschensis and R. rickeri: In the species tree, these form a supported clade that is sister to the stigmatica group, and nested within Rhyacophila. In the Bayesian inference R. rickeri is sister to all other Rhyacophila, while H. kuldschensis remains nested within Rhyacophila Additionally, R. vibox was not placed with support in the species tree (bs = 0.39).

The results from the Bayesian multistate analysis showed that acceptance rates of proposed rate changes stabilized around 38% indicating sufficient mixing of the MCMC chain. The resulting probabilities for ancestral character states in the binary coding (no gills/gills) were the same as those for the sum of all individual gill states vs no gills. In accordance with our sampling design and hypotheses we focused on the ancestral character state estimations in the R. vulgaris group (sensuSchmid 1970; Fig. 3, Clade D). The R. vulgaris group comprises four different gill types, making it the most diverse in terms of gill types known to date. However, we also estimated ancestral character states along the phylogeny backbone leading to the vulgaris group. The analyses indicate a gillless ancestor for the Rhyacophilidae and the genus Rhyacophila with Himalopsyche nested within (Fig. 3, nodes I–IV). The ancestor for Himalopsyche cannot be estimated based on the single representative of the genus. Nodes V–VII also show gillless larvae as the most likely ancestral state. This pattern changes in clade C, when the estimate of ancestral states along the backbone first shifts to four digitate gills at nodes VIII–XI and then to single-tufted gills at nodes XII–XIII (Fig. 3). Comb-like gills are estimated to arise with node XIIa. Based on ancestral state estimations, the evolution of gills shapes for the vulgaris group would thus be no gills > four digitate gills > single-tufted gills. Under these estimates the simple single digitate gills would have arisen either from a gillless ancestor or an ancestor with four digitate gills (Fig. 3, nodes VIII–IX), whereas the complex comb-like gills either arose from an ancestor with four digitate gills or single-tufted gills (Fig. 3, nodes XI–XII, XIV).

Of the 10 recognized gill types in Rhyacophila it appears that only single digit gills (R. loxias and R. laevis) evolved more than once. While tufted gills also evolved multiple times, they did so with different configurations in the vulgaris group (single tuft), brunnea group (three tufts) and grandis group (two tufts) (Giersch & Wissemann 2012). Ancestral character estimations indicate convergent evolution in both cases.

The distribution of river zonation occurrences for individual gill types is presented in Figure 4. The data used are summarized in Supplementary Material C. Species with single-tufted gills have the broadest distribution and extend furthest downstream. Their distribution is significantly different from that of species without gills or species with single digitate gills (ANOVA, F_44,4 = 5.175; p = 0.00167; Tukey HSD: R. sensu stricto vs Hyporhyacophila p = 0.0034; R. sensu stricto vs Prosrhyacophila p = 0.042). While there appears to be a trend of species with increasing structural complexity and gill surface area to occur further downstream, differences between other gill types are not significant.