A previous survey in several populations of Uruguay reported the presence of nine different morphs similar to Hyalella curvispina. In autumn 2020, we visited nine country points where these different morphs came from and collected 8 of them. Sampling localities were: San José (H8), Colonia (H1), Durazno (H3), Colonia Rosell y Rius (H2), Lavalleja (H4), Paso de los Toros (H5), Batoví (H7), Achar (H6), and two localities of Montevideo: Facultad de Ciencias (FC) and Montevideo type locality for Hyalella curvispina (MVD) (Fig. 1 and Table S1). We collected the specimens with a water net with a diameter of 20 cm and a mesh width of 250 µm. All samples were stored in liquid N₂ at collection sites and preserved there until the molecular procedure. Two adult males from each site were dissected and morphologically evaluated under a stereoscope instead of glass to confirm the presence of the characteristics of each morph according to the collection site.

Figure 1. 

Map of collection sites of the different Hyalella morphs in Uruguay. H1: Colonia (34°26.0501′S 57°49.4717′W); H2: Colonia Rossell (33°10.9503′S 55°44.363′W); H3: Durazno (33°24.0584′S 56°31.18′W); H4: Lavalleja (34°30.4393′S 55°22.3598′W); H5: Paso de los Toros (32°45.4339′S 56°31.8643′W); H6: Achar (32°23.8234′S 56°9.4333′W); H7: Batoví (31°52.9283′S 56°0.7117′W); H8: San José (34°18.65′S 56°52.7′W); FC: Montevideo (34°52.8334′S 56°7.05′W); MVD (34°50.3666′S 56°16.05′W): Montevideo (Hyalella curvispina).

We obtained total genomic DNA extractions from five animals for each of the nine morphs previously collected, following standard protocol for precipitation of proteins with salts and DNA with ethanol (modified from Miller et al. 1988). We used only the pereion and pleon of the animal to avoid contamination of possible microorganisms adhering to the mouthparts, the concentration of DNA extractions were estimated in the nanodrop spectrophotometer (NanoDrop, Thermo Scientific, USA). We analyzed the genes COI, 12S, 28S and H3 and amplified them by PCR using the primers LCO and 5587F (Stutz et al. 2010), 12Samphi-f and 12Samphi-r (Rodrigues 2016), Rnest and Fnest (Stutz et al. 2010) and H3af and H3ar (Colgan et al. 1998) respectively. The amplification were carried out in a PX0.2 Thermal Cycler (Thermo Electron Corporation) in a total volume of 22 μl containing 10µl of kit GoTaq Hot Start Green Master Mix (Promega), 2 µl of each primer (10 mM) and 8µl DNA dilution (1/100). We included negative controls in all cases, replacing the DNA dilution with water. Sequences of all morphs were obtained, except for COI. In particular, we got an average of 2 COI sequences for five morphs and H. curvispina, an average of 3 sequences of 12S, and one sequence of H3 and 28S for all morphs and H. curvispina of the two Montevideo localities (FC and MVD).

We amplified a 369 base pair (bp) fragment of the mitochondrial cytochrome c oxidase subunit I (COI) gene was using the primers LCO and 5587F (Stutz et al. 2010). The PCR program consisted of initial denaturation for 2 min at 95 ° C followed by 36 cycles of 30 s at 95ºC, 30 s at 45ºC, 1 min at 72ºC, and a 7 min extension step at 72ºC (modified from Stutz et al. 2010). We amplified the 12S mitochondrial ribosomal gene using the primers 12Samphi-f and 12Samphi-r (Rodrigues 2016) to obtain a partial sequence of 452 bp. The PCR cycle was an initial denaturation for 5 min at 96°C, followed by ten cycles of 30 s at 96°C, 60 s at 55°C (decreasing annealing temperature by 1°C/cycle during seven cycles and then decreasing to 47°C in the last cycle), and 60 s at 72°C, followed again by 30 cycles of 30 s at 96°C, 60 s at 45°C and 60 s at 72°C, with a final extension of 5 min at 72°C (modified from Rodrigues 2016). We amplified the ribosomal nuclear gene 28S with primers Rnest and Fnest (Stutz et al. 2010) to obtain a partial sequence of 605 bp. The PCR program consisted of a 1 min denaturation step at 94ºC, 39 cycles of 1 min at 94ºC, 1 min at 51°C, 1 min at 72°C, and a 5 min extension step at 72°C (Stutz et al. 2010). Finally, we amplified a 332 bp of the H3 protein-coding gene H3af and H3ar (Colgan et al. 1998). The PCR program consisted of initial denaturation for 5 min at 96°C, followed by 30 cycles of 30 sec at 96°C, 45 s at 50°C, and 60 s at 72°C, with a final extension of 5 min at 72°C (Rodrigues 2016). PCR products were checked at electrophoresis with agarose gel (0.8% in TBE 1X), stained with Ethidium Bromide, and visualized with an ultraviolet light source. PCR products of the expected size without secondary bands were purified and automatically sequenced (Sanger method, Macrogen Inc., http://www.macrogen.com) from both ends.

We edited obtained sequences manually using the PROSEQ 3.2 program (Filatov 2002), considering that each change corresponded to well-defined peaks in the chromatogram. In the case of protein-coding genes (COI and H3), we checked the absence of stop codons or indels that could modify the reading frame to ensure no pseudogene was present. In nuclear genes, we annotated the polymorphisms with an IUPAC ambiguity code.

We compared the obtained sequences with sequences from different species of the same genus reported in Genbank (see accession numbers and other details in Supplementary Table S2). All genes were aligned independently using the MUSCLE algorithm with the MEGA X program (Kumar et al. 2018). In all cases, a sequences from Platorchestia sp. (Table S2), a genus closely related to Hyalella, was used as the outgroup taxon.

For each gene, nucleotide frequencies, variable sites, and parsimony-informative (Pi) sites were estimated using MEGA X software. We calculated global nucleotide distances and pairwise distances between morphs and species with the K2P substitution model (Witt et al. 2006). The most suitable nucleotide substitution model for each gene was assessed with JMODELTEST v2.2.10 (Darriba et al. 2012), with Akaike Information Criterion (AIC), and with Bayesian Information Criterion (BIC) (Posada and Crandall 1998). Additionally, we estimated phylogenetic reconstruction with MEGAX by applying different criteria with 1000 bootstrap pseudoreplicates and pairwise deletion options (i.e., eliminate all positions with less than 95% site coverage). We carried out the reconstruction by maximum likelihood with the substitution model suggested by JMODELTEST. The initial tree(s) for the heuristic search were automatically obtained by applying Neighbor-Joining and BioNJ algorithms to a matrix of pairwise distances estimated using the Maximum Composite Likelihood distance. For maximum parsimony, we used the Subtree-Pruning-Regrafting algorithm with search level 1, in which the initial trees were obtained by the random addition of sequences (10 replicates). For neighbor-joining algorithms, evolutionary distances were computed using the Maximum Composite Likelihood method and are in the units of the number of base substitutions per site.

Bayesian phylogenetic methods were also performed in BEAST v.2.6.3 (Suchard et al. 2018) and carried the option *Beast, to estimate the species tree by considering the information of all markers. For each gene, we used the best-fit nucleotide substitution model estimated by JMODELTEST v2.2.10 (Darriba et al. 2012). We explored the nucleotide substitution saturation for each molecular marker using the Xia test in DAMBE v. 7.2.144 (Xia 2018). For COI, we separated the analysis and eliminated the third position of the codon. We partitioned the COI gene, analyzing the first and second codon positions separated from the third due to saturation of this marker (see discussion below). We linked the mitochondrial genes, the two partitions of COI and 12S, kept the other two genes independent, and selected the Yule process. It used the MCMC length of 100.000.000 generations, sampling every 10.000 with a burn-in of 20% of the trees sampled. Using TRACER v. 1.6 (Rambaut 2018), we assessed the resulting log files and corroborated that the Effective Sampling Size (ESS) values were higher than 200 to ensure adequate sampling and convergence. We created a maximum clade credibility tree in TREE ANNOTATOR v2.6.3 (Suchard et al. 2018), and this tree was visualized and edited in FIGTREE v. 1.4.4. (Rambaut 2018). Due to the difficulties in obtaining the COI gene for all the Uruguayan samples, and because the evolutionary history of mitochondrial genome is also in 12S (by linkage, as the mitochondrial genome does not recombine), we estimate two species trees, with and without COI, using the vouchers sequenced for four and three genes respectively. Finally, we applied a maximum likelihood approximation implemented in IQ-TREE v1.6.12 (Nguyen et al. 2015), with a partition by gene and unliked edges, and an automatic selection of the best substitution model for each gene, and 1000 pseudoreplicates of ultrafast bootstrap (Minh et al. 2013).

We applied three species delimitation methods for each gene independently and for all of them concatenated (with and without COI). We used the ABGD method developed by Puillandre et al. (2012) on the webserver http://wwwabi.snv.jussieu.fr/public/abgd/abgdweb.html. After sequence alignment, we computed a matrix of pairwise distances using the K2P model (Kimura 1980). We used Pmin = 0.001 and Pmax = 0.1 and X = 0.5. We also applied a recently developed method, still based on pairwise genetic distances, but whose implementation provides a score for each defined partition and overcomes the challenge of a priori defining p, the ASAP (Puillandre et al. 2021). ASAP was run in webserver https://bioinfo.mnhn.fr/abi/public/asap/asapweb.html, with K2P model and default options. Additionally, we used PTP and bPTP on the webserver https://species.h-its.org/ptp. PTP is a model for delimiting species on a rooted phylogenetic tree developed by Zhang et al. (2013). It models speciation or branching events in terms of the number of substitutions, so it only requires a phylogenetic input tree. And bPTP is an updated version of the original maximum likelihood PTP that adds Bayesian support values to delimited species on the input tree. Higher BS value on a node indicates all descendants from this node are more likely to be from one species.

Because we faced difficulties in amplifying and sequencing COI sequences, we recovered 11 COI sequences that include only Hyalella curvispina from FC locality in Montevideo, and five morphs. We could amplify and sequence 30 sequences of the 12S gene, ten 28S and H3 of H. curvispina, and all morphs (Table 1).

We summarized the sequence length obtained from each marker and the number of parsimony-informative sites in Table 2. Although 28S was the longest sequence (605 bp), the proportion of informative sites was low (approx. 18%). We analyzed shorter sequences from two mitochondrial markers, 369 and 452 bp, COI and 12S, respectively. They presented a high proportion of informative sites from the point of view of parsimony (39 and 43%, respectively). Histone H3 was the smallest fragment (332bp) and was the one that presented the highest proportion of parsimoniously informative sites, although, in absolute terms, there were only a few (22 sites) (Table 2).

Most of the morphs have different alleles/haplotypes, but there were shared haplotypes in all cases, generally for geographically close localities. In all cases: i) the H. curvispina samples (topotypes) are identical to those of FC locality close to 20 km (Montevideo), except for the H3 gene, ii) localities H1 and H8 with a distance of 120 km between them (Colonia and San José, respectively) are identical except for the COI, iii) localities H5 and H6 with a distance of 67 km between both (Paso de los Toros and Achar) except for COI. For the COI, one of the samples from the locality H1(1) (Colonia) shares a haplotype with H. curvispina (type locality and another). However, for mitochondrial 12S, the H1 and H8 morphs (from Colonia and San José, respectively) are identical and different from those of H. curvispina (both localities). In contrast, the morphs H5 and H6 of the localities (Paso de los Toros and Achar) are identical.

The 28S marker has the lowest variation, and the samples from populations H5, H6 (Paso de los Toros and Achar), and H. curvispina (two locations) on the one hand, as well as those from H1, H8, and H3 (Colonia, San José, and Durazno), are identical.

For histone H3 the samples from H1, H3, H4, H5, H6, H7 and H8 (Colonia, Durazno, Lavalleja, Paso de los Toros, Achar, Batoví and San José) are identical, while Hyalella curvispina from the two localities of Montevideo do not share haplotype. We kept identical variants in the analyzes to estimate the species tree with all of them (Table 3).

The substitution models selected by the AIC criterion were HKY+G for the COI and 12S, the TVM+G for the 28S, and K80+G for the H3. The substitution models selected by BIC were TIM3+G for COI and 12S and TMV+I+G for 28S and TPM1uf+G for H3.

Pairwise genetic distances are summarized in Table 4 and Supplementary Table S3, S4, S5, S6. A low genetic distance is observed for all markers (average 10, 4, 1 and 1%, for COI, 12S, 28S, and H3, respectively) among the Uruguayan morphs, with specimens H4(1) and H4(2v) (Table S3) being from Lavalleja, the much more divergent morph concerning the rest of the Uruguayan samples (19, 8, 2 and 1% for COI, 12S, 28S, and H3, respectively). As expected, for all markers except COI, the average genetic distances between Hyalella spp. and Platorchestia sp. used as outgroup taxon are much greater than the distances within the genus Hyalella (Table 4). For the COI, the average genetic distance between all the species evaluated and H. azteca (the most basal species of all those analyzed according to the study of Zapelloni et al. (2021) (at the genomic scale) is practically the same as the distance to the outgroup taxon (24 and 25%, respectively).

The average distances among Uruguayan samples are moderate (0.10, 0.04, 0.01, 0.01 for COI, 12S, 28S, and H3, respectively). We found the highest genetic distance between the H4 (specimens H4(1) and H4(2v)) and the rest (0.19, 0.08, 0.02 and 0.01, for COI, 12S, 28S, and H3, respectively), similar to the interspecific distances between the Uruguayan samples (excluding specimens H4(1) and H4(2v)) and other species such as H. kochi 2319B.

For all markers, the sequences showed little substitution saturation. For each marker Iss is significantly lower than Iss.c for both symmetric and asymmetric topologies (these were always lower than the firsts). For COI sequences, for the three codon positions taken together (Iss: 0.22 < Iss.c: 0.76 (symmetric topology) Iss.c: 0.52 (asymmetric topology), degree of freedom (DF): 114, p < 0.0001), or for the 1st and 2nd position of the codon taken together (Iss: 0.05 < Iss.c: 0.90 (symmetric topology) Iss.c: 0.76 (asymmetric topology), DF: 75, p < 0.0001). For the 12S mitochondrial gene sequences (Iss: 0.13 < Iss.c: 0.68 (symmetric topology) Iss.c: 0.36 (asymmetric topology), DF: 256, p < 0.0001). For 28S sequences (Iss: 0.07 < Iss.c: 0.72 (symmetric topology) Iss.c: 0.41 (asymmetric topology), DF: 521, p < 0.0001) and for histone H3 sequences (Iss: 0.04 < Iss.c: 0.53 (symmetric topology) Iss.c: 0.41 (asymmetric topology), DF: 95, p < 0.0001).

No gene analyzed independently showed solvency in the resolution of the most nodes (see phylogenetic reconstructions (Figs 2–5, and in Supplementary Figs S1–S8). On the contrary, most of the nodes have low to moderate support values. Only some were high and varied with the reconstruction method used. However, except for the H3 gene, the COI and 12S mitochondrial genes and the 28S nuclear gene show that the Uruguayan samples are paraphyletic. There is a monophyletic group formed by Uruguay, H. montforti 2015 2D, and H. kochi (4747, 2319B, 3TK27), which includes two groups: one wider group that includes H. curvispina and has a higher affinity with H. montforti 2015 2D and H. kochi 4747, 2319B, 3TK27, and the other with fewer samples, sister to the previous one. In these tree topologies, we observed that the samples annotated as H. kochi are polyphyletic. The H3 gene presents a very low variation (showed in other studies) and could not resolve any node with bootstrap values from moderate to high. We do not recover the monophyly of the Uruguayan samples, nor the close relationship between these and H. montforti 2015 2D or some samples of H. kochi 4747, 2319B, 3TK27.

Figure 2. 

Phylogeny of Hyalella reconstructed by Maximum Likelihood using all Uruguayan specimens sequenced for COI (369 bp), 28 Hyalella sequences from North and South America, and an outgroup taxon (Platorchestia japonica). Uruguayan samples, all in “curvispina complex,” are denoted in grey color. Bootstrap values are next to the nodes.

Figure 3. 

Phylogeny of Hyalella reconstructed by Maximum Likelihood using all Uruguayan specimens sequenced for 12S (452 bp), 21 Hyalella sequences from North and South America, and an outgroup taxon (Platorchestia parapacifica). Uruguayan samples, all in “curvispina complex,” are denoted in grey color. Bootstrap values are next to the nodes.

Figure 4. 

Phylogeny of Hyalella reconstructed by Maximum Likelihood using all Uruguayan specimens sequenced for 28S (605 bp), 27 Hyalella sequences from North and South America, and an outgroup taxon (Platorchestia japonica). Uruguayan samples, all in “curvispina complex,” are denoted in grey color. Bootstrap values are next to the nodes.

Figure 5. 

Phylogeny of Hyalella reconstructed by Maximum Likelihood using all Uruguayan specimens sequenced for H3 (332 bp), 18 Hyalella sequences from North and South America, and an outgroup taxon (Platorchestia pacifica). Uruguayan samples, all in “curvispina complex,” are denoted in grey color. Bootstrap values are next to the nodes.

The Bayesian species tree considering all markers, each with its most appropriate substitution model, resolves most nodes with higher supports of the posterior probability. Among them, we can highlight 1) the monophyly of the Uruguayan samples together with two samples of H. kochi 4747, 2319B, and H. montforti 2015 2D, the basal position of H4(1) concerning this clade, and the reciprocal monophyly of the Uruguayan samples without H4(1) on one side and of H. montforti 2015 2D with two samples of H. kochi (2319B and 4747) on the other, 2) the basal position of H. azteca and H. armata 26-2A within the species of the genus analyzed, 3) the monophyly of two other groups of species H. nefrens 2310E + H. hirsuta 30-5C + H. neveulemairei 30-5D + H. kochi AP18, and H. cajasi EC3-1 + H. tiwanaku 2304 (Fig. 6), although with very low statistical support (0.39).

The species tree obtained from the two nuclear genes and only one mitochondrial (only 12S and excluding COI) has the same general topology and higher posterior probability values. The information from the Uruguayan sample H7(1) (not recovered with COI) is incorporated and show higher affinity with two samples of H. kochi (4747, 2319B) and H. montforti (2015 2D) than with the majority group of Uruguayan samples (Fig. 7).

Figure 6. 

Phylogeny of Hyalella reconstructed by Bayesian analysis. Samples of H. curvispina FC1, H1(1), H2(1), H3(1), H5(1) and H. sp.1 H4(1) are Uruguayan samples. Eleven Hyalella sequences from North and South America and outgroup taxon (sequences of different specimens of the genus Platorchestia sp.) were included. The consensus tree is based on 1758 bp from COI, 12S, 28S, and H3 concatenated datasets. Posterior probabilities are noted next the nodes. Clades from A to F were defined in Zapelloni et al. (2021); Clade G is proposed in the present study.

Figure 7. 

Phylogeny of Hyalella reconstructed by Bayesian analysis. Samples H. curvispina FC1, MVD, H1(1), H2(1), H3(1), H5(1), H6(1), H7(1), H8(1) and H. sp.1 H4(1) are Uruguayan samples. Eleven Hyalella sequences from North and South America and outgroup taxon (sequences of different specimens of the genus Platorchestia sp.) were included. Shown is the consensus tree based on 1389 bp from 12S, 28S and H3 concatenated datasets. Posterior probabilities are noted next to the nodes. Clades from A to F were defined in Zapelloni et al. (2021); Clade G is proposed in the present study.

Maximum likelihood reconstructions yielded similar results (Fig. S9 and Fig. S10) for four or three genes, respectively). The major difference in both cases is the absence of reciprocal monophyly between groups H. nefrens 2310E + H. hirsuta 30-5C + H. neveulemairei 30-5D + H. kochi AP18, and H. cajasi EC3-1 + H. tiwanaku 2304.

The results of molecular species delimitation are shown in the Supplementary Fig. S11. The results for the different genes analyzed independently are similar, quite different, and COI was the gene that yielded the most different results between methods. Among the four methods tested, ASAP found more groups in all genes, and differed most from the other methods.

In general, H. armata 26-2A, H. azteca, H. cajasi EC3-1, H. tiwanaku 2304, and H4(1), a sample of the “curvispina complex”, stand out as well-differentiated species with high statistical support. For all genes and methods, sample H4(1) was always clearly different from the rest of the Uruguayan samples, with the exception of the ABGD method applied to COI and concatenated with COI. Besides, a group of different species is considered by these methods as a single species: H. kochi AP18, H. neveulemairei 30-5D, H. nefrens 2310E, and H. hirsuta 30-5C. In this group H. kochi AP18 is differentiated analyzing H3 and 28S (by ABGD) and H. hirsuta 30-5C analyzing COI (by PTP and bPTP). ASAP excludes from this group H. kochi AP18 and H. hirsuta 30-5C analyzing COI and 12S. The Uruguayan samples, except H4(1), generally cluster as a single species in COI, 12S, and 28S with ABGD, PTP and bPTP, and the concatenate without COI (by PTP and bPTP).

The phylogeny obtained incorporates representatives of the genus Hyalella from Uruguay, including H. curvispina, into the general phylogeny previously proposed for the genus. As expected, the general topology concerning the other species coincides with that previously reported by Zapelloni et al. (2021), which used identical Hyalella sequences. However, we used only those vouchers sequenced for all genes. We recovered all proposed clades by Adamowicz et al. (2018) (except B, which we did not include) as monophyletic with high support. In addition to previous studies, H. azteca and H. armata 26-2A (Clade C) are more basal than the southern clades. The closer relationship between these species presents a posterior probability value of less than 0.5 so we cannot trust the phylogenetic relationships inferred at these nodes, as Zapelloni et al. (2021) show. The South American species are clustered together and include different clades and Uruguayan samples. However, the relationships between groups A, E, D, and F could not be resolved with significant statistical support (except for basal group C).

The “curvispina complex” forms a monophyletic group with H. montforti 2015 2D and H. kochi 4747, 2319B, 3TK27, both corresponding to clade E from the northern Altiplano (Zapelloni et al. 2021). However, H. kochi is also elsewhere in the phylogenetic tree, so some of these vouchers considered as H. kochi are probably cryptic species. Then, the close association between H. curvispina and H. kochi (sensu stricto) should be taken with caution. In any case, all the samples associated with Uruguayan Hyalella come from Peru, which is very distant geographically. This study is the first one that includes Hyalella species to the east of the Andes, i.e., “curvispina complex”. The close link between the samples from Peru and Uruguay may reflect the scarce knowledge about the genus. It would be necessary to include many more South American representatives (in addition to those already included from Peru, Bolivia, and Chile). In the future, this phylogenetic information could shed light on the historical biogeography of Hyalella and the palaeoclimatic history of the continent.

On the other hand, we found high genetic variations between H. curvispina and associated morphs and at least one provisional new species. We consider that specimens collected at the type locality of H. curvispina, which also have the expected morphology, are indeed topotypes of this species. Other samples, collected in distant lakes (except H7(1) and H4(1) + H4(2v)), show higher affinity with that species, with high statistical support in the species phylogeny that includes the COI (96%). Indeed, many of these morphs together with H. curvispina are genetically identical in some of their markers. The genetic differentiation found among most of them was moderate, in the range expected for intraspecific differentiation, and consistent among most genes (Supplementary Tables S3, S4, S5, S6). Thus, we propose that all these morphs are part of the H. curvispina variation generated in the region. All the Hyalella specimens found in Uruguay have morphological characteristics that define the “curvispina complex”, they are: smooth body surface, presence of curved setae on inner side of inner ramus of uropod I and sternal gills present in segments 2 to 7. On the other hand, we found one morph, H4 specimens H4(1), H4(2v) and H4(3) (only for 12S gene), which show greater differentiation than the rest of the morphs (Fig. 2, 3, 4, 5 and Supplementary Tables S3, S4, S5, S6). Although the COI showed significant differentiation between the different morphs, this is not consistent with the information provided by another mitochondrial marker, 12S, nor by nuclear markers (see below). We propose maintaining the complex name Hyalella curvispina for most of these forms, and suggest one new provisional species corresponding to specimens H4(1), H4(2v) and H4(3)) (Lavalleja locality) named H. sp.1 in phylogenetic trees. The morphological characteristics that differentiate H. sp.1 from H. curvispina are subtle, namely, telson with three strong setae distributed in two groups in the apical margin and inner median face of propodus in gnathopod 1 with an oblique row of 9–10 short pectinate setae.

In addition, all Uruguayan samples are more closely related to clade E (97%), so they could be considered part of clade E. Genetic distances between clades measured as K2P for 28S are in the range of 1.1% to 6.4% (Adamowicz et al. 2018). The H7(1) sample belongs to clade E, with a genetic distance to it of 0.4%. The remaining samples from Uruguay, except for H. sp.1, belong to this clade with a distance of 1.1%. In addition, we suggest a new clade, clade G, formed by H. sp.1; the genetic distance between this new clade and clade E is 1.8%.

The information provided by each of the markers independently, and the markers as a whole, suggests a geographical differentiation within the H. curvispina clade (part of group E). We now considered “H. curvispina clade”, the clade including the type locality and related localities, excluding H. sp.1 and the sample H7(1), with the closest populations being more genetically associated. However, the relationship between these “geographic groups” is still uncertain because of low statistical support. And shallow supports would reflect a recent differentiation. A multilocus approach, including thousands of markers (e.g., obtained from RADseq and NGS approximations), will probably be helpful to resolve the critical fine-scale aspect of phylogeography needed to ascertain further details of this differentiation.

In the locality of Lavalleja, we collected three samples, two of them with a higher genetic divergence (in mitochondrial and nuclear markers) regarding the variation of most Uruguayan variants/morphs and preliminary morphological evaluation shows differences between this sample and H. curvispina (Waller, data not shown). We suggest that this sample, i.e. specimens H4(1), H4(2v) and H4(3), would be considered a different species H. sp.1. The other one H4(2) presents low genetic divergence and morphologically corresponds to H. curvispina. Molecular species delimitation analyses strongly support this conclusion. Different methods identify the sample H4(1) as an independent species for all genes except COI and their concatenation. Since this sample belongs to a new species, it confirms sympatric species living together in the same pool. Several authors have been observed sympatric distributions in the two species complexes evaluated at the molecular level. In particular, the cryptic species of H. azteca from North America (Wellborn and Cothran 2004; Witt and Hebert 2000) and in Hyalella species from Brazil (da Silva Castiglioni and Bond-Buckup 2008a, 2008b, 2009; González et al. 2006). On the other hand, sample H7 shows greater affinity to two samples of H. kochi (4747, 2319B) and to H. montforti 2015 2D than to the rest of the Uruguayan samples. Although H7 was sequenced for mitochondrial 12S and nuclear markers and not for COI due to technical problems, we propose this sample be a new species. However, molecular species differentiation analyses do not discriminate it as a different species, but it is integrated into the E clade. It would be necessary to review the systematic of group E as a whole and establish clearer species boundaries or suggest synonymy to H. montforti (by the principle of priority of the zoological nomenclatural code). Overall, these results highlight the relevance of including molecular systematics studies in determining the genus Hyalella.

The Species Screening Threshold criterion (SST) (Witt et al. 2006) has been used in the molecular systematics of the genus for the last two decades (Dionne et al. 2011; Vergilino et al. 2012) and considers that Hyalella species diverge from each other by 3.75% for COI sequences using the K2P distance (Kimura 1980). Following this criterion, the presence of cryptic species could be considered for the H. curvispina species complex, particularly six species of the morphs (excluding H7 not sequenced for COI). However, our analysis performed with a 369 bp of COI does not provide solid phylogenetic support for this delimitation. Strikingly several recent Hyalella COI barcoding studies have used even shorter fragments, sometimes only 300 to 400 bp in length, for specimen identification (Dionne et al. 2011; Major et al. 2013; Stutz et al. 2010).

On the other hand, comparing the information provided by COI with that from various markers, both mitochondrial and nuclear, reveals the specific difficulties associated with this marker. Firstly, the genetic distances in the COI marker are practically the same between different levels of variation. In particular, values in the order of 4 to 29% are found between poorly differentiated species within the “curvispina complex”, as well as between highly differentiated species (some species of the “curvispina complex” with other Hyalella species), or between Hyalella species with the outgroup taxon This condition may be due to saturation (i.e., multiple substitutions at the same site in a sequence leads to underestimation of actually occurring mutations) and leads to homoplasy and an underestimation of divergence times between haplotypes observed, particularly for older phylogenetic events (e.g., Wilke et al. 2009). However, when analyzing the degree of saturation with the DAMBE program, low saturation levels were observed for all markers, including COI without the third codon position. The absence of significant-high saturation for COI could be due to the sequence size (369 bp), as saturation in COI for Hyalella has been recorded at sequence sizes larger than 500 bp (Major et al. 2013; Worsham et al. 2017; Zapelloni et al. 2021).

In turn, the 12S gene gives us different information to the COI, although both are mitochondrial genes and are physically linked because this genome does not recombine (Meyer 1993; Saldamando and Marquez 2012). Although both markers share the evolutionary history, the 12S gene is the most conserved within the mitochondrial genome (Arif and Khan 2009). Thus, unlike the COI, 12S might not be saturated (Major et al. 2013; Witt and Hebert 2000; Worsham et al. 2017) and offers more accurate information at this level of comparison (Arbogast et al. 2002; Major et al. 2013). On the other hand, pairwise genetic distances with all markers except COI (with nuclear genes, 12S, or even the concatenated construct of all genes) show different sharpie levels of variation within the genus (i.e., intra- and interspecific) and among genera (between Hyalella and outgroup taxon). Pairwise genetic distances calculated with COI (Table 4) do not show this pattern, and variations within and among genera are of similar magnitude. Thus, we believe that for Hyalella, 12S is more reliable than COI and that the SST criterion should be applied to another marker.

We found little differentiation in the markers assessed in this study, both nuclear and mitochondrial 12S. However, mitochondrial differentiation at the COI level is high (in fact saturated), which has also been reported in studies conducted for other complexes (Major et al. 2013; Witt and Hebert 2000; Worsham et al. 2017). These results may reflect the current differentiation processes of the genus across regions, in principle dominated by colonization and extinction events (Zapelloni et al. 2021). Similarly (Duan et al. 2000), and despite the difference in geographical scale of the studies, unique variants and high levels of variation are also found between populations 50–200 km apart. This differentiation is associated with the geographic variation. Closer populations are more phylogenetically linked, suggesting some connectivity between populations and diversification in the presence of gene flow within H. curvispina (sensu latu).

Also, in agreement with previous studies about other species complexes in the genus, we found that the “curvispina complex” is paraphyletic respect to the species H. kochi 4747, 2319B, 3TK27 (sensu latu) and H. montforti 2015 2D. These results suggest that adaptive and morphological convergence in this group is high (Adamowicz et al. 2018), and the “curvispina complex” is no exception. In particular, morphological variation does not match genetic differentiation, which may be related to the recurrent selection of similar morphologies in the face of the same ecological different challenges. In contrast to the genetic phylogeny, all samples from Uruguay, including H. sp.1 and sample H7(1), are more similar morphologically to each other and H. curvispina than to H. kochi (2319B, 4747 and 3TK27) and H. montforti 2015 2D. Habitat specialization and trophic regimes could explain the convergence of morphotypes observed in Hyalella specimens from Uruguay. Hyalella curvispina has varied feeding habits: shredder, predator, scrapper, and collector-collector (Cummins et al. 2005; Giorgi and Tiraboschi 1999; Saigo et al. 2009; Wantzen and Wagner 2006). They are also food for other macroinvertebrates, fish, amphibians, and birds (Colla and César 2019). As a defense mechanism against predation pressure, Hyalella species occupy different habitats (Wellborn 1995). They inhabit a variety of freshwater environments such as lakes, ponds, and streams, clinging to vegetation and burrowing in bottom sediments, where they are important members of the benthic fauna (da Silva Castiglioni and Bond Buckup 2008B; Grosso and Peralta 1999; Wellborn 1995).

On the other side, the two phylogenetically most closely related species to the H. curvispina complex (i.e., H. kochi and H. montforti) share some morphological characteristics with it. Regarding H. kochi both species have a smooth body surface, the inner face of propodus of gnathopod 1 with seven setae, the presence of curved setae in the inner ramus of uropod 1. However, the main characteristic that distinguishes these two species is the presence of six pairs (from segment 2 to 7) of sternal gills in H. curvispina, while H. kochi has only five pairs (from segment 3 to 7), and the consistency of this character makes it relevant in the evolutionary relationships within the genus (González and Watling 2001). Compared with H. montforti, both species have six pairs of sternal gills and curved setae in the inner ramus of uropod 1. Still, the main characteristic that distinguishes these two species is the body with dorso-posterior flanges on pereon segment 7, pleonite 1, 2, and 3 in H. montforti. Since these forms have a common ancestor and live in different environments, it is reasonable to think that the common morphological characteristics reflect phylogenetic inertia, while others would be local adaptations.