The specimen was swept from the riparian vegetation using a hand net.

A pair of legs from the holotype of R. lignumvallis sp. nov. were removed, and the tissue was incubated overnight in 60 µl TNES lysis buffer (100 mM Tris-HCl, 25 mM NaCl, 10 mM EDTA, 1% SDS) and 8 µl Proteinase K (20 mg/ml). For DNA binding and cleanup, 75 µl custom speed-bead suspension was added (Sera-Mag SpeedBeads Carboxylate, hydrophobic, Cytiva; see Rohland and Reich 2012), incubated for 15 min on a rotating shaker, and the beads were washed twice with 75% ethanol after the supernatant had been removed and discarded. The DNA was eluted in 1X TE.

DNA sequences were generated for the cytochrome c oxidase subunit I barcoding region (mtCOI, 658 bp) using primers LCO1490-L and HCO2198-L (Nelson et al. 2007). Polymerase chain reactions (PCR) were run on a Mastercycler Pro S (Eppendorf, Hamburg, Germany) in reactions containing 1X MyTaq Reaction Buffer, 0.4 µM of each forward and reverse primer, 0.5 U MyTaq DNA Polymerase, 1 µl DNA and nuclease-free water to fill up to a total volume of 10 µl. Reaction conditions were 1 min at 95 °C for initial denaturation followed by 35 cycles of 20 s at 95 °C (denaturation), 30 s at 45 °C (annealing), and 30 s at 72 °C (extension). The reaction ended with a final extension for 5 min at 72 °C. PCR products were visualized on agarose gels and purified using a modified ExoSAP protocol with Exonuclease I (20U/µl) and Fast AP Themosensitive Alkaline Phosphatase (1U/µl; both ThermoFisher Scientific, Vilnius, Lithuania). DNA sequences were generated at the Laboratory Centre of the Senckenberg Biodiversity and Climate Research Centre using a 3730XL DNA Analyzer (Applied Biosystems).

Genomic DNA was taken from the above DNA isolates, quantified using a Qubit 4.0 fluorometer with a 1x dsDNA HS Assay Kit (ThermoFisher Scientific, Waltham, USA), and sheared to a target fragment size of 350 bp using a Bioruptor Pico (Diagenode, Seraing, Belgium). Genomic libraries were prepared from 27.4 ng sheared gDNA using the NEBNext Ultra II DNA Library Preparation Kit for Illumina (New England Biolabs, Ipswich, MA, USA), following the manufacturer’s manual. Adapters were diluted 1:10 as recommended for low input libraries, and size selection was omitted due to the low DNA amount. A dual indexing PCR was run for 7 cycles on a Mastercycler (Eppendorf, Germany) using NEBNext Multiplex Oligos for Illumina (Dual Index Primers Set 1; New England Biolabs, Ipswich, MA, USA). After cleanup, library integrity was verified on a 2200 TapeStation with a High Sensitivity D1000 Tape (Agilent, Santa Clara, CA, USA), and shipped for 150 bp paired-end sequencing (ordering 30 Gbp output) on a partial lane of an Illumina NovaSeq 6000 platform (San Diego, CA) at Novogene (Cambridge, UK).

Raw reads are deposited in the National Center for Biotechnology Information’s Sequence Read Archive (NCBI SRA) under the accession number SRR22799047 under Bioproject PRJNA899095.

The final mtCOI alignment included 71 sequences: 68 sequences from 16 species and 2 morphospecies in the Rhyacophila tristis group. Rhyacophila italica Moretti, 1981, Himalopsyche kuldschensis (Ulmer, 1927), and H. triloba (Hwang, 1958) were included as outgroups. All sequences were generated for this manuscript, except R. bosnica Schmid, 1970 (MK211322), H. kuldschensis (KX143534), and H. triloba (KX295339), which were retrieved from GenBank. The barcode region is 658 bp in length. However, some of our sequences were incomplete, and their lengths were between 576 and 658 bp. Therefore, the final alignment was completed with Ns on both ends, reaching 3.69% of missing data. Sequence specimen data and GenBank accession numbers are summarized in Table 1.

The maximum likelihood tree was produced using IQ-TREE v.2.1.3 (Minh et al. 2020), using the command iqtree2 -s RhyacophilaLignumvalleMS_658bp.fasta -B 1000 -bnni -alrt 1000 --prefix RhyacophilaLignumvalleMS_658bp. The TIM2+F+I+G4 nucleotide substitution model was selected using ModelFinder (Kalyaanamoorthy et al. 2017). Statistical support for the tree topology was assessed with the ultrafast bootstrap approximation (UFboot) (Hoang et al. 2017). Clades with UFboot ≥ 95% are considered well-supported. Additionally, each bootstrap tree was optimized with a hill-climbing nearest neighbor interchange (NNI) search (flag -bnni in the command above) based on the corresponding bootstrap alignment to prevent overestimating UFboot branch support values, as recommended by Hoang et al. (2017).

The consensus tree was visualized and edited in TreeViewer v.2.2.0 (Bianchini and Sánchez-Baracaldo 2024). Additional aesthetic edits were made in Adobe Illustrator CS6.

After quality control with FastQC v.0.11.9 (Andrews 2019), raw reads were trimmed for low-quality regions, adapter sequences, and overrepresented k-mers using autotrim.pl v.0.6.1 (Waldvogel et al. 2018) and Trimmomatic v.0.39 (Bolger et al. 2014) with the adapter_all.fa of Trimmomatic and the following settings ILLUMINACLIP:2:30:10:8:true, SLIDINGWINDOW:4:20, MINLEN:50, and TOPHRED33. Unpaired reads were discarded and paired reads were checked for contamination using Kraken v.2.0.9 (Wood and Salzberg 2014).

Genome size was estimated using a method based on k-mer distribution. For this, k-mers were counted with JELLYFISH v.2.3.0 (Marçais and Kingsford 2011) using jellyfish count -C -s 1000000000 -F 2 and a k-mer length of 21 based on the raw sequence reads. A histogram of k-mer frequencies was created with jellyfish histo and used for analysis with the online web tool GenomeScope v.2.0 (Ranallo-Benavidez et al. 2020) using the following parameters: k-mer length = 21, ploidy = 2, max k-mer coverage = 10000.

The mitochondrial genomes were first assembled with the raw reads using NOVOplasty v.4.2 (Dierckxsens et al. 2016) using the following parameters: type = mito, genome range = 12000–22000, k-mer = 33, max memory = 100, read length = 150, insert size = 300, platform = illumina, paired = PE, insert size auto = yes. The partial sequence of the cytochrome c oxidase subunit I (COX1) gene of Rhyacophila fasciata Hagen, 1859 (MT559357.1) was used as seed input. All other parameters were kept as default. In addition, we used a second mitogenome assembler MitoZ v.2.3 (Meng et al. 2019). For this purpose, the raw data was subsampled to 10,000,000 reads using seqk and then used as input for MitoZ assemble with the following parameters: genetic_code 5, clade Arthropoda, fastq_read_length 150, insert_size 300, run_mode 2, filter_taxa_method 1, requiring_taxa ‘Arthropoda’. Annotation of tRNAs, rRNAs, and protein-coding genes was done for the best mitogenome assembly of each species with MitoZ v.2.3 using the module annotate with genetic_code 5 and clade Arthropoda. Both mitogenome assemblies were aligned to the complete mitogenome of R. quadrifida Sun & Yang, 1995 (OL678049.1) and R. kando Schmid, 1970 (OL678048.1) with MAFFT in Geneious Prime v.2022.1.1 (Biomatters Ltd.) to set the correct start position and manually curate the control-region. The mitochondrial genome assembly was deposited in GenBank under the accession OQ984043.

Nuclear genome assembly was conducted in Spades v.3.14.1 (Bankevich et al. 2012) with the default settings. Scaffolds smaller than 500 bp and those with blast hits to the mitochondrial genome assembly were filtered out. Assembly statistics were calculated with Quast v.5.0.2 (Gurevich et al. 2013) and completeness was assessed via screening for single-copy orthologs with BUSCO v.5.2.2 (Manni et al. 2021) using the endopterygota_odb10 dataset. As an additional quality control, trimmed reads were mapped back to the assembly with bwa-mem v.0.7.17-r1188 (Li 2013) with parameters -a -c 10000, and the back-mapping rate was calculated with qualimap v.2.2.1 (Okonechnikov et al. 2015). To check for potential contamination, taxon-annotated GC-coverage (TAGC) plots were generated with BlobTools v.1.1.1 (Laetsch and Blaxter 2017). For this purpose, the bam file resulting from the back-mapping analysis was converted to a blobtools readable cov file with blobtools map2cov. Taxonomic assignment for BlobTools was conducted with blastn v.2.10.0+ (Camacho et al. 2009) using -task megablast and -e-value 1e-25. The blobDB was created and plotted with the cov file and blast hits. The nuclear draft genome assembly was deposited in GenBank under accession JAPMAE000000000. The DNA barcode region was extracted from the genome assembly and aligned to the traditionally sequenced mtCOI sequences and showed 100% identity to the larvae included by Engelhardt et al. (2011).

Genomic methods were identical for R. tsurakiana (see Suppl. material 3).

The holotype specimen was prepared and examined following standard methods for ethanol-preserved material (Blahnik and Holzenthal 2004; Blahnik et al. 2007). Forewing length was measured from base to apex with a microscale (BioQuip Products, Rancho Dominguez, California, USA).

The abdomen was removed from the specimen, soaked in 85% lactic acid, and heated to 99 °C for 60 min to dissolve internal soft tissues. The macerated tissues were then flushed out of the abdomen with a syringe. The holotype was examined on an Olympus SZX10 stereomicroscope, and pencil sketches of the genitalia were prepared using a Leitz Dialux 20 compound microscope outfitted with a drawing tube. Pencil sketches were scanned and placed in an Adobe Illustrator CS6 document as a template for vector illustrations. Morphological terminology follows Schmid (1970) for the male genitalia, Holzenthal et al. (2007) for wing venation, and Ivanov (1990) for setal warts.

After collapsing clades with less than 70% bootstrap support, the species in the Rhyacophila tristis group were placed in a polytomy, and sister to R. italica (Fig. 2). The first clade of the polytomy includes R. pubescens, R. tsurakiana, and the new species. The holotype specimen (marked with an asterisk in Fig. 2) was included in a highly supported clade (100% bootstrap support) with three Rhyacophila larvae from Corsica. Based solely on larval identification, these were originally considered to be R. pubescens (Engelhardt 2009; Engelhardt et al. 2011). It now seems clear that these are larvae of R. lignumvallis sp. nov.

Figure 1. 

Type locality of Rhyacophila lignumvallis sp. nov. on the island of Corsica (France).

Figure 2. 

Phylogenetic relationships in the Rhyacophila tristis group based on the mtCOI barcode region and estimated by maximum likelihood, as implemented in IQ-TREE v.2.1.3 (Minh et al. 2020). Nodal support was calculated using the ultrafast bootstrap (UFboot) approximation (Hoang et al. 2017); nodes with UFboot values greater than 95% are considered well supported. Branches with support values of less than 70% were collapsed. Branch lengths are measured in nucleotide substitutions per site (see scale bar for reference). The branch leading to Rhyacophila was shortened (total branch length was 0.12 substitutions per site). Numbers in parentheses correspond to GenBank accession numbers. The holotype and the Corsican larvae are shaded.

The second clade in the polytomy includes several species from the tristis subgroup (R. orghidani Botosaneanu, 1952, R. cibinensis Botosaneanu & Marinkovic-Gospodnetic, 1967, R. margaritae Kumanski, 1998, R. bosnica, and R. obtusa Klapalek, 1894), and two species from the pubescens subgroup (R. pirinica Kumanski, 1980 and R. producta McLachlan, 1879). These specimens were collected in Albania, Austria, Bulgaria, and Romania. The third clade in the polytomy exclusively includes species from the tristis subgroup (R. aquitanica McLachlan, 1879, R. carpathica Botosaneanu, 1995, R. trescavicensis Botosaneanu, 1960, R. tristis Pictet, 1834, R. vranitzensis Botosaneanu & Marinkovic-Gospodnetic, 1967, and Rhyacophila sp.). These specimens were collected in Austria, Bosnia and Herzegovina, Bulgaria, France, Italy, Montenegro, Romania, and Spain.

Illumina sequencing resulted in 212,866,450 raw reads with a data amount of 31.9 Gbp for R. lignumvallis sp. nov. After trimming and contamination filtering, 173,132,236 reads (22.2 Gbp) were kept. The Genomescope2 analysis revealed a genome size of 699,853,381 bp and heterozygosity of 20% (see Suppl. material 2).

The NOVOplasty mitogenome assembly resulted in three contigs (18,087 bp, 1,404 bp, 238 bp) that could not be circularized. Therefore, the 15,623 bp long contig obtained by MitoZ was chosen for annotation. The annotation of the mitogenome revealed all expected 13 protein-coding genes and both rRNAs and 23 tRNAs.

The nuclear genome assembly of R. lignumvallis sp. nov. contains 206,802 scaffolds with a total length of 644 Mb, an N50 of 5.6 kb, and a GC of 30%. The BUSCO search with 2,124 Endopterygota orthologs resulted in 82.5% BUSCOs; of these, 47.9% were complete (47.4% single, 0.5% duplicated), and 34.6% were fragmented. 96.3% of the reads were mapped back to the original assembly. Blobtools detected no contamination in the assembly for R. lignumvallis sp. nov. (see Suppl. material 2). However, some contamination was detected by NCBI using the improved FCS-GX screen according to https://github.com/ncbi/fcs (see Suppl. material 2).

The genomic characterization of R. tsurakiana is included in Suppl. material 3.

The authors have declared that no competing interests exist.

No ethical statement was reported.

Financial support from the LOEWE Centre for Translational Biodiversity Genomics, Frankfurt am Main, Germany (LOEWE-TBG).

ER: investigation, data curation, formal analysis, writing – original draft, writing – review and editing, visualization. WG: investigation, writing – original draft, writing – review and editing, visualization. JH: data curation, formal analysis, writing – original draft, writing – review and editing. JS: investigation, writing – original draft, writing – review and editing. SUP: study design, data curation, investigation, writing – original draft, writing – review and editing.

Ernesto Rázuri-Gonzales https://orcid.org/0000-0002-7554-0816

Wolfram Graf https://orcid.org/0000-0001-6559-0644

Jacqueline Heckenhauer https://orcid.org/0000-0001-8771-9154

Julio V. Schneider https://orcid.org/0000-0002-9823-6569

Steffen U. Pauls https://orcid.org/0000-0002-6451-3425

All of the data that support the findings of this study are available in the main text or Supplementary Information.