Baetis rifensis sp. nov. is the first representative of the Baetis fuscatus group to be described from the Maghreb. It was collected from streams in the Rif region of northern Morocco. All species of the B. fuscatus group are morphologically very similar, with slight differences in colour. Thus, in addition to morphological description, species delimitation based on genetic evidence was carried out. The new species was compared with other members of the B. fuscatus group from the Palaearctic region.

COI, Maghreb, mayflies, morphology, Rif, taxonomy

Baetidae is the most species-rich family of mayflies, with approximately 1,100 species in 114 genera (Sartori and Brittain 2015; Jacobus et al. 2019; Kaltenbach et al. 2021), representing roughly one-third of all mayfly species globally. The family is cosmopolitan, except for New Zealand (Gattolliat and Nieto 2009; Gattolliat et al. 2023). Taxonomic and hydrobiological studies are constantly adding to the list of baetid species, but knowledge remains incomplete. The mayfly fauna in Morocco comprises 54 species, with half of them belonging to Baetidae. A comprehensive checklist of these species has been recently published; however, some species are still awaiting description (El Alami et al. 2022a; Kaltenbach et al. 2022; Gattolliat et al. 2023).

Thanks to the arduous efforts of specialists in Morocco, the country has been extensively investigated since 1990. The Rif, the solitary mountain range that emerged from the alpine orogeny and situated in the northernmost region of Morocco, has received particularly comprehensive research attention. Encompassing an area of 30,000 km² (Fig. 1), the Rif stretches approximately 90 km from north to south and spans 340 km from east to west (El Alami and Dakki 1998; Errochdi et al. 2012; El Yaagoubi et al. 2023).

Figure 1. 

Distribution of Baetis rifensis sp. nov. (yellow star: type locality).

The genus Baetis Leach, 1815 presents a complicated taxonomic history, as most Baetidae with hind wings and double intercalaries in the fore wings were historically assigned to this genus. Some of the taxa were treated in the last 30 years as species groups or subgenera of Baetis, including taxa that are no longer considered to belong to the genus Baetis, such as Labiobaetis Novikova & Kluge, 1987 (Monaghan et al. 2005; Webb 2013; Kluge and Novikova 2016). The old concept of Baetis is now considered to be polyphyletic, as indicated by its division into multiple genera (e.g. Waltz and McCafferty 1987; Webb et al. 2018), which are partly based on species groups proposed for European species (Müller-Liebenau 1969). Overall, Baetis fuscatus and Baetis vernus groups are included in the genus Baetis s.s. according to the classification widely used in Europe (Bauernfeind and Soldán 2012; Webb et al. 2018).

The Baetis fuscatus group was defined by Müller-Liebenau (1969); it encompasses the European species B. beskidensis Sowa, 1972, B. fuscatus (Linnaeus, 1761), and B. scambus Eaton, 1870 (Bauernfeind and Soldán 2012), and the North American species B. flavistriga McDunnough, 1921, B. intercalaris McDunnough, 1921, B. notos Allen & Murvosh, 1987, B. phoebus McDunnough, 1923, and B. rusticans McDunnough, 1925 (Webb et al. 2018). Morphological differences between the species are very limited. The general colouration of the imagos as well as the pronotum pattern of the larvae are often used for separating the species. The B. fuscatus group was reported in the Eastern Palaearctic region, notably in locations such as Mongolia (Braasch 1983; Landa and Soldán 1983), Korea (Bae and Soldán 1997; Bae and Yoon 1997), China (Quan et al. 2002), and Siberia (Novikova and Kluge 1997). Unfortunately, no molecular data are currently available for these materials. Additionally, confidently attributing these populations to a specific species remains a challenging task. While a comprehensive revision of the Eastern Palaearctic species is highly recommended, it should be noted that such an endeavour falls beyond the scope of the present study.

In North Africa, the B. fuscatus group was first discovered in Algeria, where it was identified as B. bioculatus. It was later also reported from Morocco generally under B. fuscatus or B. gr. fuscatus (Thomas 1998; El Alami et al. 2022a). Baetis fuscatus group was mentioned from the Rif (El Alami 1989, 2002; El Bazi et al. 2017; Khadri et al. 2017; Guellaf et al. 2021), while in contrast, it seems to be absent in the Central Plateau and Oriental Morocco (El Agbani et al. 1992; Mabrouki et al. 2017; El Alami et al. 2022a).

In this study, we describe a new species, Baetis rifensis sp. nov., based on larvae collected from Rif streams. We used an integrative approach combining morphological and genetic evidence to separate the Rif populations from other Western Palearctic species.

By integrating all data at our disposal, we provide the distribution of Baetis rifensis sp. nov. in the Rif area since 1990. Specimens were collected from aquatic habitats (intermittent, ephemeral, or perennial streams) using a kick sampler during field excursions. Larvae were captured in different seasons. The data for this study came from fieldwork conducted by the first author from 2019 to 2023 and from previously collected material by El Alami between 1989 and 2018 from nearly identical sampling points.

All samples were preserved in 70% ethanol solution. The larvae were dissected and mounted on microscopic slides for morphological study under an Olympus SZX7 stereomicroscope in the Laboratoire Écologie, Systématique, et Conservation de la Biodiversité, Université Abdelmalek Essaâdi (LESCB).

Photographs of captured larvae were taken using a Canon EOS 6D camera and processed using Helicon Focus v. 5.3 (Helicon Soft Ltd.) and Adobe Photoshop Lightroom v. 5 (Adobe Systems Inc.). Adobe Photoshop Elements v.13 was used to enhance each image. Photographs of body parts of the larvae were taken with an Olympus BX43 microscope equipped with an Olympus SC50 camera and processed with Olympus software Cell Sense v. 1.3. The distribution map was generated using ArcGIS software (ESRI, Inc). The GPS coordinates of the sample locations are given in the sections of examined type materials.

To complement morphological investigations, we sequenced a 658-bp fragment of the mitochondrial gene cytochrome oxidase subunit 1 (COI) for five specimens of B. rifensis sp. nov., following the non-destructive DNA extraction procedure described in Vuataz et al. (2011). The polymerase chain reaction (PCR), purification, and sequencing steps were conducted according to the methodology described in El Alami et al. (2022b). Forward and reverse sequencing reads were assembled and edited in Codon Code Aligner v. 10.0.2 (Codon-Code Corp.) and aligned using MAFFT (Katoh et al. 2019) with default settings as implemented in Jalview v. 2.11.2.6 (Waterhouse et al. 2009). The number of parsimony-informative sites of the alignment was calculated in MegaX (Kumar et al. 2018; Stecher et al. 2020). To augment our dataset, we initially downloaded all COI sequences associated with members of the B. fuscatus group available on the GenBank database as of 11 May 2023, resulting in a total of 385 records. Additionally, we obtained the sequences accessible on the BOLD Systems data portal as of the same date and selectively retained only those that were not shared with GenBank, yielding an additional set of 74 sequences. We also included three additional sequences (two from Spain and one from Portugal) obtained from the unpublished European mayfly project FREDIE (https://wp.fredie.eu/). To reduce computational demand and improve gene tree readability, we then manually excluded GenBank/BOLD sequences obtained from specimens collected outside the Western Palearctic region. This selection process was conducted after confirming that the COI sequences of B. rifensis sp. nov. were clearly distinct from the removed sequences (data not shown). A total of 64 sequences remained for further analyses, comprising five newly generated sequences (Table 1), 55 sequences from GenBank (seven from Sroka 2012; two from Cardoni et al. 2015; two from Corse et al. 2017; 15 from Morinière et al. 2017; seven from Múrria et al. 2017; seven from Tenchini et al. 2018; 11 from Behrens-Chapuis et al. 2021; four unpublished iBOL data releases), one from BOLD (GMGMQ2692-18), and three from the project FREDIE (ES029_SR10F10; ES030_SR10G01; PT002_SR18E12).

To explore and visualize the COI evolutionary divergence, we employed pairwise genetic distances and gene-tree approaches. COI pairwise distances were calculated using the dist.dna function from the ape 5.7-1 package (Paradis and Schliep 2019) in R v. 4.2.3 (R Core Team 2023), selecting the raw model and the pairwise.deletion option, corresponding to uncorrected p-distances (see Srivathsan and Meier 2012) with missing data removed in a pairwise way. Mean, minimum, and maximum distances within and between putative COI species, referred to as molecular operational taxonomic units (MOTUs) hereafter, were calculated using the ddply function from the plyr 1.8.8 package (Wickham 2011). The assignment of COI sequences to MOTUs was determined based on the results of the ASAP species-delimitation analysis (as described below). Prior to reconstructing the COI gene tree, the best evolutionary model (GTR+Γ+I) was selected based on the second-order Akaike information criterion (AICc; Hurvich and Tsai 1989) implemented in JModelTest v. 2.1.10 (Darriba et al. 2012) with five substitution schemes, six gamma categories, and default values for other parameters. To account for different substitution rates among COI codon positions, we analyzed our data set in two partitions, one with first and second codon positions, and the other with third positions (1 + 2, 3). A Bayesian-inference analysis was performed using BEAST v. 1.10.4 (Suchard et al. 2018) on the CIPRES Science Gateway 3.3 (Miller et al. 2010). The input BEAST file was generated in BEAUTi (Suchard et al. 2018), incorporating the selected evolutionary model and partition scheme described above. A relaxed molecular clock model (uncorrelated lognormal) and a UPGMA starting tree were used, with default settings for other parameters. Two independent Markov-chain Monte Carlo (MCMC) chains were run for 30 million generations, logging parameters every 1000 generations. Convergence of the MCMC runs was visually verified in Tracer v. 1.7.2 (Rambaut et al. 2018). The log and tree files from the independent runs were combined using LogCombiner v. 1.10.4 (Suchard et al. 2018), after discarding the initial 10% of trees as burn-in, ensuring that all parameters reached effective sample size values > 200. The maximum clade-credibility tree was obtained using TreeAnnotator v. 1.10.4 (Suchard et al. 2018) with default settings. Visualization and editing of the tree were conducted in iTOL v. 6.5.7 (Letunic and Bork 2021).

Finally, we applied two contrasting single-locus species delimitation methods to our COI dataset: the distance-based ASAP (Assemble Species by Automatic Partitioning; Puillandre et al. 2020) and the tree-based mPTP (multi-rate Poisson Tree Processes; Kapli et al. 2017) approaches. ASAP, an improved version of the ABGD (Automatic Barcode Gap Discovery; Puillandre et al. 2012) approach, was employed using the ASAP webserver (https://bioinfo.mnhn.fr/abi/public/asap/asapweb.html) to estimate the most probable number of MOTUs based on our COI alignment. We calculated genetic distances using simple p-distances and selected the species delimitation hypothesis associated to the highest barcode gap width (W) among the three partitions sharing the same best asap-score. The mPTP method, which is a multi-rate extension of the PTP (Poisson Tree Processes; Zhang et al. 2013), takes advantage of phylogenetic differences within and between species based on the number of substitutions obtained from a phylogenetic tree. We conducted mPTP using the web service available at https://mptp.h-its.org, using the BEAST COI gene trees as input (see above).

The material is deposited in the collections of the Laboratoire Écologie, Systématique, et Conservation de la Biodiversité (LESCB) in Tétouan, and the Muséum cantonal des sciences naturelles in Lausanne (MZL).

The widespread adoption of integrative approaches for studying Moroccan mayflies has the potential to enhance our understanding of their species diversity in the coming years. Past dubious identifications could be revised, and the description of new endemic species could be facilitated, as in the recently described species of Prosopistoma maroccanum (El Alami et al. 2022b) and Centroptilum alamiae (Kaltenbach et al. 2022). Furthermore, environmental factors go a step further in explaining the detected pattern of current mayfly distribution, as expected from ecosystems with high spatial and temporal heterogeneity, such as Mediterranean rivers. This prompts the need to address and fill the existing geographic and taxonomic gaps by directing future sampling missions to the most under-prospected locations across the entire country. This will help to remove any uncertainties in the actual identification of the genus Baetis.

We express our gratitude to Céline Stoffel for her invaluable assistance in molecular biology and the preparation of COI barcodes. We are very grateful for the insightful comments and helpful recommendations offered by the reviewers, Luke M. Jacobus and Lucas Lima.