Morphological and molecular species boundaries in the Hyalella species flock of Lake Titicaca (Crustacea: Amphipoda)

The Hyalella species diversity in the high-altitude water bodies of the Andean Altiplano is addressed us-ing mitochondrial cox1 sequences and implementing different molecular species delimitation criteria. We have recorded the presence of five major genetic lineages in the Altiplano, of which one seems to be exclusive to Lake Titicaca and nearby areas, whereas the rest occur also in other regions of South America. Eleven out of 36 South American entities diagnosed by molecular delimitation criteria in our study are likely endemic to the Titicaca and neighbouring water bodies. We have detected a remarkable disagreement between morphology and genetic data in the Titicacan Hyalella , with occurrence of several cases of the same morpho-species corresponding to several Molecular Operational Taxonomic Units (MOTUs), some even distantly related


Introduction
Ancient lakes -those with an uninterrupted history of more than 100,000 years (Gorthner, 1994) -are remarkable by their exceptional species richness compared to other continental water bodies.In such type of lakes, some groups of organisms are prone to occur in the form of large "swarms" of closely related endemic species (Brooks, 1950;Fryer, 1991).For instance, Lake Baikal harbours at least 257 species and 74 subspecies of amphipod crustaceans, all but one being endemic to the lake, representing around 14% of the known global diversity of freshwater amphipods (Kamaltynov, 1999).Even more remarkable are the ca.1,000 and 600 endemic species of cichlid fish estimated to live in the East African rift valley lakes Nyasa and Victoria, respectively (Fryer, 2000).How this extraordinary diversity has originated and diversified in a limited space poses fundamental questions on the mechanisms of speciation, community structure and coexistence of closely related species (Cristescu et al., 2010).
A remarkable aspect of these swarms of lacustrine species is the frequent occurrence of a mismatch between species delimited based on DNA sequences and those demarcated following more traditional, morphology-based criteria (Kroll et al., 2012;Naumenko et al., 2017).Such taxonomic decoupling has been recorded to occur also in some of the large radiations of oceanic island terrestrial invertebrates and may be related to the fast speciation rate that concurs on such insular environments (Monaghan et al., 2006;Fryer, 1996;Cristescu et al., 2010).Factors such as genetic introgression, retention of ancient polymorphisms, morphological convergence among distant genetic lineages through selection, or through a combination of selection and hybridization, could be possible causes of the failure of both classical morphological and DNA barcoding approaches to resolve the taxonomy of many of these rapidly radiating lineages (Monaghan et al., 2006).This issue poses a challenge for the species identification and inventorying of the rich biota present in ancient lakes.
Lake Titicaca, at 3,806 m of altitude, is the only ancient lake present in South America.Although its age has not been precisely determined, it is assumed to date back to between 3 and 2 Million years (Ma) (Kroll et al., 2012).In any event, it is not older than the formation of the northern Andean Altiplano where it is located, which has been dated at 5.4 ± 1.0 Ma (Kar et al., 2016).The Titicaca is placed in a vast endorheic area (∼200,000 km2) which harbours also Lake Poopó and the salt flats of Coipasa and Uyuni (fig.1A, B).Like many other ancient lakes, the Titicaca has a complex palaeo-environmental history (Cristescu et al., 2010).Thus, the different Altiplano basins conformed a single hydrological unit from the Early to Middle Pleistocene (Baker & Fritz, 2015;Nunnery et al., 2019), when water level experienced at least three major expansions (high stands) of up to 114, 104 and 54 m above present lake level (Fornari et al., 2001).Furthermore, the lake has passed through several periods of retraction coincident with global interglacial periods, when its level FIGURE 1 A, map of South America showing location of the Andean Altiplano; B, approximate area of the Altiplano showing placement of its main water bodies; C, sampling sites placed outside Lake Titicaca; D, sampling sites at Lake Titicaca itself.Numbers identify sampling stations.See supplementary tables S1 and S2 for precise information on sampling sites.Maps were produced using R ggmap (Kahle & Wickham, 2013) and Google Maps (Google Inc., Mountain View, CA).
Downloaded from Brill.com02/23/2021 11:52:37AM via IMEDEA Instituto Mediterraneo de Estudios Avanz dropped well below its outlet threshold, leading to a closed-basin configuration and an increase in water salinity (Fritz et al., 2012).In the central and southern Altiplano -where Lake Poopó and the salt flats of Coipasa and Uyuni are located -, there is also evidence of the existence of palaeo-lakes at least since the Late Pleistocene.These lakes originated by waters overflowing from the Titicaca basin after the progressive erosion of the threshold between the northern and central Altiplano, a process that was not completed until ca.60,000 yr BP (Fritz et al., 2004).In summary, the Altiplano lakes have a complex history of periodical basin connections, disconnections, episodes of desiccation, and abrupt changes in water salinity and depth, all of which may have exerted an important effect on the pattern of species diversification in the area (Nunnery et al., 2019;Fornari et al., 2001).
The molecular phylogeny of the Hyalella of Lake Titicaca has been recently studied by Adamowicz et al. (2018) using mitochondrial cox1 and nuclear ribosomal 28S sequences.They concluded that they are polyphyletic, identifying at least five major distinct evolutionary lineages within the lake, all of them with an age estimated to be much older than the lake itself (20-12 Ma vs. maximum age of the lake 5.4 Ma).Only two of these clades seem to have further diversified by an apparently recent intra-lacustrine diversification (Adamowicz et al., 2018).
In the present study, we bring together DNA-sequence data and morphological observations derived from new material collected at Lake Titicaca and its main satellite lakes, plus at many other high-altitude water bodies of the Altiplano (see fig. 1) to explore the Hyalella species diversity and their diversification pattern.Expanding the sampling scheme to cover as many localities as possible around the Titicaca is crucial to elucidate the modes of diversification of its amphipod species flock since many of these high-altitude lakes were captured by the Titicaca in former periods of high-water level (see above).We make use of an extensive cox1 dataset resulting from merging our sequences to the ones available in BOLD (www.boldsystems.org)to explore the Hyalella species diversity in the Altiplano using mPTP, GMYC and ABGD species delimitation criteria.The species delimitation proposals are contrasted with the morphological variation displayed by Titicacan Hyalella, exploring the following questions: (1) have the armoured bodies of some Hyalella species Downloaded from Brill.com02/23/2021 11:52:37AM via IMEDEA Instituto Mediterraneo de Estudios Avanz evolved only once or did they arise independently in several Titicacan lineages?(2) how congruent are the different molecular delimitation methods implemented and how do they relate to the morphospecies known to occur in the Andean Altiplano?Our study also aims to build a genetic reference cox1 library as a platform to be used in future morphological and evolutionary studies of this amphipod lineage endemic to an ecologically complex, ancient lake system.

Samples
We collected Hyalella specimens at 51 different sites in Lake Titicaca and at 35 sites from 30 additional water bodies of the Altiplano and nearby areas of Perú and Bolivia (see fig. 1C, D, and supplementary tables S1 and S2 for the location of sampling sites and composition of the Hyalella assemblage per site).Furthermore, we used also three samples of H. cajasi Alonso & Jaume, 2017, from southern Ecuador (collected at high-altitude lakes of El Cajas Massif; Azuay), and one of H. franciscae González & Watling, 2003, from southernmost Chile (Madre de Dios Island; Magallanes Region), all of them out of the study area.Specimens were collected directly with hand-held plankton nets or with a small dredge thrown either from the shore of the lake or from a boat, and preserved in 96% ethanol.Individuals were first sorted out and identified based on morphological characters and classified as "morphospecies".To explore the correspondence between morphospecies and Molecular Operational Taxonomic Units (MOTUs) we used ten taxa easily and unambiguously diagnosable based on morphology (see below for a key to all Hyalella species from Lake Titicaca and other water bodies of the Altiplano area).Unfortunately, we failed to collect any representative of the species triplet H. echinus/ crawfordi/ gauthieri, easily distinguishable by the common display of combined midline and lateral armature on body tergites (see Coleman & González, 2006).The species H. tiwanaku, H. kochi and H. cuprea are likely to include hidden diversity and were used only as a preliminary morphospecies reference.The species considered were as follows: H. longipes (eleven spines on dorsal midline; fig.2K); H. neveulemairei (six spines on dorsal midline; fig.2F); H. longipalma (five spines on dorsal midline, one on each pereionites 6-7 and pleonites 1-3, respectively; fig.2B); H. solida (also five spines on dorsal midline, but this time one on each pereionites 5-7 and pleonites 1-2, respectively; fig.2A); H. armata (one spine laterally-directed on each pereiopodal coxal plates 1-4; fig.2G, H); H. knickerbockeri (two spines on dorsal midline, one on each pleonites 1-2, respectively; fig.2I); Hyalella n. sp. 1 (body smooth, compact, with shortened antennules and antennae; fig.2J); Hyalella n. sp. 2 (body smooth, with medial margin of merus and carpus of pereiopods 3-4 fringed with long setae); H. nefrens (three spines on dorsal midline, one on each pereionite 7 and pleonites 1-2, respectively, and mandible incisor smooth); H. montforti (four dorsal spines, one on each pereionite 7 and pleonites 1-3, respectively, and with male first uropod with a modified spine on endopod; fig.2C); and H. latimanus (also with four dorsal spines as H. montforti, but with male first uropod lacking modified spine on endopod).

Molecular procedures
Genomic DNA was individually purified from whole specimens using the Qiagen DNeasy Blood & Tissue kit (Qiagen, Hilden, Germany) following the manufacturer protocol.Voucher specimens and DNA aliquots are deposited at the Instituto Mediterráneo de Estudios Avanzados (IMEDEA, Spain).A total of  197 specimens were sequenced for a fragment of the mitochondrial cytochrome c oxidase subunit 1 gene (cox1) using the primer-pair LCO/ HCO (Folmer et al., 1994), and a subset consisting of 69 individuals for a fragment of the nuclear gene Histone H3 using the primer pair H3aF/H3aR (Colgan et al., 1998) (GenBank accession numbers MN582154 -MN582350 and MN582351 -MN582419, respectively; see supplementary table S3).The cox1 sequences obtained were distributed as follows: 183 individuals from the Titicaca and nearby areas in Perú and Bolivia, 13 from Ecuador and one from Chile (see fig. 1C, D, and supplementary tables S1 and S2 for details on sampling sites).For phylogenetic and species delimitation analyses our sequence data set was merged with the cox1 sequences of South American Hyalella available at the BOLD database (http://www.boldsystems.org;project "TTKK"; included partially in the study of Adamowicz et al., 2018).The resulting dataset represents an extensive South American Hyalella cox1 repository of 1,139 sequences from Perú and Bolivia (mostly collected at Lake Titicaca and areas nearby), 122 from Argentina, 65 from Chile and 13 from Ecuador.Sanger sequencing followed standard protocols described elsewhere (Bauzà-Ribot et al., 2011).Sequences were filtered using the Perl script uniqHaplo.pl(Takebayashi, 2015), retaining unique haplotypes for subsequent analyses.The resulting dataset consisted of 560 Hyalella cox1 sequences and was used for subsequent phylogenetic and species delimitation analyses.The histone H3A sequences rendered only nine haplotypes with a total of 23 parsimony-informative nucleotide positions and a remarkably low phylogenetic signal (see results).

Phylogenetic analyses
The 560 cox1 mitochondrial DNA sequences were trimmed to start and end with complete triplets, aligned at the amino acid level using MAFFT 7 (Katoh & Standley, 2013) in SeaView (Gouy et al., 2009) to inspect for the absence of indels and stop codons, and finally backtranslated to nucleotides preserving the alignment by codons, resulting in final nucleotide alignment of 648 bp.Histone H3 sequence alignment was also performed with MAFFT7 using the default FFT-NS-1 algorithm.As no robust support for a reciprocal monophyletic relationship of North and South American Hyalella was obtained in preliminary analyses, seven Platorchestia spp.cox1 sequences retrieved from GenBank (accession numbers KC578469, KC578473, KC578478, KC578488 and KC578494, KC578495, MG936436) were used as outgroup in the mitochondrial phylogenetic analyses.Cox1 phylogenetic trees were built under a maximum likelihood (ML) framework in IQ-TREE (Nguyen et al., 2014) and under Bayesian analyses with BEAST v1.8.4 (Drummond & Rambaut, 2007).The best partition schemes and models were explored in PartitionFinder (Lanfear et al., 2012).The optimal partitioning strategy and evolutionary models consisted of subdividing the cox1 data set by codon positions with the models TRN+I+G, GTR+G and GTR+I+G, respectively.For Bayesian analyses, best partition schemes, two different molecular clock (strict and lognormal-relaxed) and two diversification models (constant size coalescent and Yule speciation process) were compared using Bayes Factors based on marginal likelihood estimations through path sampling and steppingstone calculations in BEAST (Baele et al., 2012).The optimal scheme consisted in partitioning the data by codon sites, and the specification of a lognormal-relaxed clock with constant size coalescence tree model (supplementary table S4).Analyses with the optimal scheme consisted of four independent runs of 200 million generations each.MCMC chain convergence and the burnin fraction (10%) was assessed in Tracer v1.

Molecular species delimitation
Species were first delimited based on the ingroup 560 cox1 unique haplotypes using genetic distances in Automatic Barcode Gap Discovery (ABGD; Puillandre et al., 2012).This approach finds recursively the slope above a cut-off value that splits the data in intra-and interspecific distances starting from an arbitrary minimum (p) and maximum (P) values.We used the default value p = 0.001, enabling the existence of nearly identical haplotypes and P = 0.1 after the maximum intraspecific distance previously found in a densely sampled amphipod species (Bauzà-Ribot et al., 2011).Distances were calculated using the K2P nucleotide substitution model with a lower cut-off value (X = 1.0).These values were used to estimate the minimum optimal gap between intra-and interspecific distances.Species boundaries and their support were also estimated with the Multi-rate Poisson Tree Processes method (mPTP), a method that fits the branching events of each delimited species to a distinct exponential distribution and does not require of any similarity threshold as input (Kapli et al., 2017).This analysis was performed based on the obtained rooted ML IQ-TREE phylogenetic tree after removing the outgroups and using two independent MCMC chains at 100,000,000 generations each with a sampling frequency of 1,000.
Finally, the Generalized Mixed Yule Coalescent model (GMYC; Pons et al., 2006;Fujisawa & Barraclough, 2013) was also used to delimitate species.This method applies a likelihood ratio test to compare a single coalescent branching rate across a clock-like tree as a null hypothesis with a model including both coalescent and Yule branching treebifurcating patterns.The R library SPLITS v1.0-19 (Ezard et al., 2009) was used for this analysis, specifying a single threshold on the ultrametric tree estimated with BEAST v1.8.4 as described above.The GMYC support value of each node representing a transition from species to population was calculated as the sum of Akaike weights of candidate delimitation models where the node is included (Fujisawa & Barraclough, 2013).A node support value of 1 means that all models tested point to the occurrence of a speciation event at that node, while lower values indicate that fewer support for it.We set an arbitrary minimum value of 0.9 to consider a GMYC entity as supported by most coalescent models.

Mitochondrial phylogeny
Both ML and Bayesian cox1 trees showed consistent topologies and yielded seven monophyletic clades with high bootstrap and posterior probability support values (called A to G from now onwards), although their inter-clade phylogenetic relationships are not well resolved with this mitochondrial marker (fig. 3 and supplementary figs S1 and S2).The Histone H3 tree topology obtained from a reduced number of the sampled individuals was less informative than the mitochondrial one but pointed to the presence of the same major clades (supple mentary fig S3).Five of these clades include taxa present in the Andean Altiplano (Adamowicz et al., 2018)

Species molecular delimitation and geographic distribution
The ABGD method identified a total of 30 MOTUs in the pool of South American Hyalella.The application of the mPTP method Downloaded from Brill.com02/23/2021 11:52:37AM via IMEDEA Instituto Mediterraneo de Estudios Avanz delimited 34 molecular entities as possible species, a number that scaled up to 56 using the GMYC algorithm, applying the rigorous option of considering only entities supported by more than 90% of the tested models.We decided to follow a conservative species delimitation scheme based on (1) consensus of the three species delimitation methods, (2) posterior probability support values of each putative species, and (3) reciprocal monophyly in the phylogenetic tree.This reduced the number of molecular entities to 36 (figs.3 and 4; supplementary fig S1) that showed a mean intergroup K2P distance of 0.171 (SD = 0.041) and an intragroup (MOTU) distance of 0.011 (SD = 0.008) (supplementary table S6).
Eight out of the 36 MOTUs identified in the analyses belong to clade A, of which all but two (MOTUs A7 and A8, limited to occur in Argentina) are circumscribed to the Titicaca and/or the Altiplano area, although some members of MOTU A1, broadly distributed across the Altiplano, occur also in the nearby Atacama region of Chile (fig. 3

Incongruence between morphological and molecular species boundaries
The mitochondrial lineages of Hyalella recognized in the Altiplano show largely overlapping geographical ranges and a considerable lack of congruence with morphological species delineations (fig. 1 and table 1).We assigned the individuals collected in the area to at least twelve already described and two previously unknown morphospecies (table 1).Two of these taxa show a broad distribution across the Altiplano and share a smooth, nonprocessiferous body, namely H. kochi and H. tiwanaku; each one of these species might represent a species complex and will not be considered further in the analyses (table 1).
Our results also show that (1) all clades but B include species with a processiferous body (figs.1).In summary, eleven out of the South American mitochondrial entities diagnosed by molecular methods in the present study appear as endemic to the Titicaca + Altiplano area (MOTUs A2-A6, C1, D1-D2, E2-E4) whereas another four seem to have expanded their range to outside the boundaries of the Altiplano (A1, B1, C2, E1), in particular to the Atacama Desert in Chile.Clades A and E are the most diversified at this high-altitude territory, including six and four genetically differentiated MOTUs, respectively.

Discussion
Our cox1 sequences coalesced in most cases with those of Adamowicz et al. (2018) data set, except for the mitochondrially delimited putative species E3, E4 and D2 from the Peruvian Altiplano lagunas Orurillo, Súchez, and Langui-Layo, respectively (stations 69, 81 and 79 in fig.1C).The connection between populations in the Altiplano and the Atacama Desert in Chile put forward repeatedly in our analysis might suggest a possible role of shorebirds in the dispersal of Hyalella species across the area (Rosine, 1956;Daborn, 1976;Swanson, 1984).Adamowicz et al. (2018) used Barcode Index Numbers (BINs), a method that implements a species threshold value of 2%, to detect 48 BINs among their South American Hyalella data set, twelve of them occurring at the Titicaca area including six uniquely sampled in the lake itself.However, studies in marine amphipods using different barcode gap Downloaded from Brill.com02/23/2021 11:52:37AM via IMEDEA Instituto Mediterraneo de Estudios Avanz detection techniques have concluded that a single threshold may be unsuitable to delimitate species in this order of crustaceans and could even be misleading (Costa et al., 2007(Costa et al., , 2009;;Tempestini et al., 2018;Jażdżewska & Mamos, 2019).Furthermore, the use of a single threshold for species delimitation can be particularly problematic in island and islandlike rapid radiating lineages (Monaghan et al., 2006).In our case, the results of the three different species delineation methods -two of them based on gene trees and one on distance matrices-differed in the number of species diagnosed in the South American Hyalella sample, with GMYC arising as the method delimiting the highest number of molecular entities.The ABGD approach produced the lowest estimate of diversity rendering a total of 30 molecular entities.The Generalized Mixed Yule Coalescent model has been shown to be useful in detecting incipient speciation events, where both intra-and interspecific distances are expected to be low, but it is known to frequently overestimate species numbers (Hendrich et al., 2010, Talavera et al., 2013, Kekkonen et al., 2015).It has been suggested that assessing the uncertainties in tree topology and the support obtained by different coalescent models in GMYC can lead to a more robust estimation of species numbers (Talavera et al., 2013, Fujisawa & Barraclough, 2013).Indeed, our Hyalella mitochondrial tree shows relatively long inter-clade but shallow within-clade branches, which is consistent with a scenario of recent diversification and the idea of "ecological opportunity" (Wagner et al., 2012).This applies in particular to clades A, D and E, that are the ones that more likely derive from incipient intra-lacustrine diversification (fig.3

and supplementary fig S1).
GMYC clearly over-splits the number of entities present in clade A, where twenty-one MOTUs were identified by this method, while ABGD and mPTP diagnosed only three and eight, respectively (supplementary fig S1).Therefore, it is preferable to apply different methodologies and to check for their congruence (Carstens et al., 2013).The Multi-rate Poisson Tree Processes method (mPTP) may be particularly well suited to our Hyalella dataset since it assumes the sampled species show different levels of intraspecific genetic diversity, arising either from differences in their evolutionary history or the uneven sampling of each taxon (Kapli et al., 2017).This method predicted the occurrence of 34 molecular entities in our sampling.We applied a conservative approach (see methods) to conclude that at least 36 molecular entities can be diagnosed in the South American Hyalella dataset, eleven of them occurring exclusively in the Andean Altiplano.
An assessment of the Hyalella species diversity present in our samples based on morphological criteria pointed to the presence of a minimum of twelve known described species and two presumably new (table 1 and  supplementary tables S1 and S2).Our results clearly show that mitochondrial sequence phylogenetic relationships are largely at odds with suggested species delineations based on morphological features, as a single MOTU, or even haplotype, is often shared by a heterogeneous array of morphospecies (see table 1).For instance, MOTU C1 is shared by H. armata (a species characterised by the display of a long transverse spine on each pereiopodal coxal plates 1-4; fig.2G, H) and Hyalella n. sp. 1 (a species devoid of any body armature and characterised by its compact body and shortened antennules and antennae; fig.2J 1).In other instances (e.g., H. neveulemairei; H. longipalma; H. latimanus) the same morphospecies includes MO-TUs falling in several of the five main clades present in the Altiplano.The widespread discordance between molecular and morphological divergence suggests an evolutionary scenario of rapid morphological and ecological differentiation accompanied by phenotypic convergence within the different lineages, a pattern common to other ancient lake radiations such as those of the cichlid fishes of the East African lakes (e.g., Meyer et al., 1990;Kocher, 2004) and the amphipod assemblage of Lake Baikal (Macdonald et al., 2005;Naumenko et al., 2017).Although the age and geological history of Lake Titicaca and Lake Baikal are notably different -as are the level of taxonomic diversity and the ancestral colonizing lineages in their respective amphipod assemblages -, a common pattern arises.It includes the remarkable level of morphological disparity attained in both amphipod assemblages; the episodic development of spiny morphologies; the presence of morphologically well-differentiated species that are hardly differentiated genetically (e.g., Parapallasea spp.and Oxyacanthus spp. at Lake Baikal) and the occurrence of cases of cryptic differentiation (e.g., Eulimnogammarus cruentus at Lake Baikal; Naumenko et al., 2017).
The phylogeny obtained shows that armoured spiny morphologies have evolved independently multiple times in the Altiplano, since they appear in all but one of the clades.Spiny morphologies are frequent among members of typical amphipod marine families such as the Iphimediidae Boeck, 1871, Dexaminidae Leach, 1814 and the Atylidae Lilljeborg, 1865, but occurs rarely among epigean continental water forms.Only some members of the Lake Baikal species-flock (Takhteev, 2000(Takhteev, , 2019)), several Caspian gammaroids (Sars, 1894(Sars, , 1896)), and the monotypic genus Fuxiana Sket, 2000 from Lake Fuxian (another ancient lake, from Yunnan, China; Sket, 2000) share with the Hyalella from the Titicaca a comparable integumentary ornamentation.Spiny morphologies have evolved independently in amphipods of Lake Baikal at least in four occasions, with at least two cases recorded of reversal to non-spiny forms (Naumenko et al., 2017).Interestingly, spinose and heavily calcified shells are also frequent among the endemic gastropods of ancient lake Tanganika, where these morphologies have appeared convergently several times.As in the case of the amphipod species flocks present in some ancient lakes, spinose shells occur frequently in marine gastropods but not in freshwater ones, a fact that has been interpreted as defensive adaptations against predators in a coevolutionary arms race scenario (Vermeij and Covich, 1978;West and Cohen, 1996).In the Titicaca, the development of this armature might be related to the predation pressure exerted by the cyprinodontid killifish endemic to the lake (Orestias Valenciennes, 1839; up to 24 spp.; Lauzanne, 1992;González & Coleman, 2002), similarly as it has been demonstrated for the amphipods Gammarus pulex and G. roeselii (Bollache et al., 2005;Copilas-Ciocianu et al., 2020).In the latter case, an extreme variation in both the number and shape of spines has been reported despite the little inter-populational cox1 differentiation, a scenario attributed to the combination of fish predation with either phenotypic plasticity or local adaptation (Copilas-Ciocianu et al., 2020).
Both our results and those of Adamowicz et al. (2018) clearly show that Lake Titicaca has been colonized by multiple South American Hyalella lineages, raising the possibility that hybridization between some of them Downloaded from Brill.com02/23/2021 11:52:37AM via IMEDEA Instituto Mediterraneo de Estudios Avanz resulted in an increase in diversification rates as has been demonstrated in the cichlid fishes of several African lakes (see below).Adamowicz et al. (2018) found a general concordance between mitochondrial cox1 and ribosomal rRNA gene (28S) phylogenies in the Titicacan Hyalella, with only one case of presumed occurrence of unidirectional hybridization.Unfortunately, our limited Histone 3 nuclear gene sequence dataset rendered a low phylogenetic signal, impeding to address a detailed analysis of the congruence between nuclear and mitochondrial gene trees.Cytonuclear discrepancies and complex relations with morphological variation are frequent issues among the cichlid fishes of the East African lakes Malawi, Tanganyika and Mweru (Egger et al., 2007;Nevado et al., 2009;Koblmüller et al., 2017;Meier et al., 2019).The origin of such discrepancies has been related to the occurrence of secondary contacts linked to lake level fluctuations among the differentiated lineages, followed by introgression or to the occurrence of incomplete lineage sorting (Koblmüller et al., 2017).
Lake Titicaca and other water bodies of the Andean Altiplano have experienced drastic changes in hydrological conditions during the geological past.Thus, retraction or even extinction of lacustrine populations might have been triggered by the increase in water salinity associated to the establishment of low water stands in the Titicaca basin.These occurred as recently as between 8,000-3,600 yr BP, when the water level receded 90 m below its present level, and especially during a period that ended up at 90,000 yr BP, when it dropped 240 m.During the latter period, the lake volume decreased by 80% and its extension and depth attained only 2,000 km2 and 40 m, respectively (D'Agostino et al., 2002).For a closed-basin interval of 9,500 years just before 2,000 yr BP, Cross et al. (2001) estimated that the maximum salinity of the lake should have increased to about half the seawater value.Although some Hyalella have been reported to thrive in waters up to 37‰ in salinity in some Canadian athalasohaline lakes (Hammer et al., 1990), the genus is mostly found at water salinities below 8 ‰ (Colburn, 1988;Galat et al., 1988;Alcocer et al., 1998).In the Altiplano area, Hyalella has not been recorded at water salinities above 2.4 ‰ (pers.obs.; Dejoux, 1993).During these past stressing periods, Hyalella might have found temporal refuge in those sectors of the Titicaca where salinity remained low due to river or groundwater discharge, or even in the assemblage of lakes and rivers of the Altiplano area, which offer plenty of opportunities to population subdivision and recolonization, enabling also the hybridization between different lineages.These episodes may have played a major role in the diversification of the Hyalella across this vast high-altitude area.
The present study examined Hyalella mitochondrial sequences distributed across Lake Titicaca and other water bodies of the Andean Altiplano.We confirmed previous results of Adamowicz et al. (2018) on the occurrence of five major monophyletic clades there, and delimited conservatively fifteen mitochondrial molecular entities within this assemblage, although they appear mostly in conflict with the morphospecies recorded in the area, which are not resolved as monophyletic groups in our cladogram.The mitochondrial genetic diversity present in the studied sequences reveals Altiplano amphipods have a complex evolutionary history that needs robust nuclear and mitochondrial phylogenetic hypotheses to explore in deeper detail which factors may have contributed to their evolutionary history.This may be achieved via the analysis of high-resolution genomic data such as complete mitochondrial genomes, orthologous nuclear genes, and Single Nucleotide Polymorphisms.
Downloaded from Brill.com02/23/2021 11:52:37AM via IMEDEA Instituto Mediterraneo de Estudios Avanz (Alas, Balaguères; France) and Edmundo Moreno (Universidad Nacional del Altiplano, Puno; Perú).Miguel Alonso kindly assisted with fieldwork in the Altiplano.We also thank the editor and two anonymous reviewers for their helpful suggestions.This work was supported by the Spanish MINECO Grant CGL2016-76164-P, financed by the Agencia Española de Investigación (AEI) and the European Regional Development Fund (FEDER).

FIGURE 3
FIGURE 3 Bayesian phylogeny and molecular species delimitation of South American Hyalella based on currently available cox1 mitochondrial data (present work and Barcode of Life DATA Systems (BOLD) repository project TTKK).Nodes with maximum nodal support are remarked with circles.Purple dots on branch tips indicate haplotypes exclusively sampled in South America outside the Altiplano area (see supplementary fig S1 for details).See main text for details on molecular species delimitation methods and results.
3 and 4); (2) occurrence of incomplete lineage sorting of mitochondrial haplotypes of morphospecies is widespread across the phylogeny; and (3) many cases of well-defined species based on morphology are in conflict with mitochondrial DNA lineages (see table1

FIGURE 4
FIGURE 4 Multidimensional scaling plot based on Kimura 2-parameters genetic distances of the 560 unique cox1 haplotypes detected in Hyalella.Dot colours refer to the main Hyalella lineages inferred in the molecular species delimitation analyses.Dots have been numbered according to their MOTU assignment.
IMEDEA (CSIC-UIB), Mediterranean Institute for Advanced Studies,C/ Miquel Marquès 21, 07190-Esporles, Balearic Islands, Spain Downloaded from Brill.com02/23/2021 11:52:37AM via IMEDEA Instituto Mediterraneo de Estudios Avanz . Members of clade D occur only at Lake Titicaca and surrounding riverine and lagoon habitats of the Altiplano, whereas clades A, B, C and E display a much broader South American distribution, encompassing Argentina and Chile aside the Altiplano of

TABLE 1
Distribution of the Hyalella morpho-species of the Titicaca/Altiplano across the different MOTUs delimited in the analyses ; see also key to species below); MOTU D1 is present in all the species with 5, 6 or 11 dorsal spines, aside of several species devoid of integumentary armature.Only the nominal species H. longipes, H. solida, H. nefrens and Hyalella n. sp. 2 display cox1 sequences falling in a single clade and MOTU (although all Downloaded from Brill.com02/23/2021 11:52:37AM via IMEDEA Instituto Mediterraneo de Estudios Avanz except the latter species share their respective MOTU with other morphospecies; see Results and table