Morphologic and genetic variation within a relict Andean catfish, Hatcheria macraei , and its relationship with Trichomycterus areolatus and Bullockia maldonadoi (Siluriformes: Trichomycteridae)

Abstract The South American siluriform fishes are found primarily in the Neotropical region, north and east of the Colorado River of Argentina, with a few relict species distributed southward and westward on both sides of the Andes Mountains. Three of these, the closely related trichomycterids Hatcheria macraei, Trichomycterus areolatus and Bullockia maldonadoi, have been subject to historical taxonomic and nomenclatural arrangements. Here, we amplify a 652-bp fragment of COI mtDNA from 55 H. macraei individuals and use publicly available Cytb mtDNA sequences of the three taxa to assess their relationship, genetic variation and haplotype distribution in relation to hydrographic basins. In addition, we extend a recent morphometric study on H. macraei by analyzing body shape in 447 individuals collected from 24 populations across their entire cis-Andean distribution. We identified some lineages previously assigned to T. areolatus that show a closer relationship to either B. maldonadoi or H. macraei, revealing new boundaries to their currently known trans-Andean distribution. We found a great morphologic variation among H. macraei populations and a high genetic variation in H. macraei, T. areolatus and B. maldonadoi associated with river basins. We highlight further integrative studies are needed to enhance our knowledge of the southern Andean trichomycterid diversity.


INTRODUCTION
The catfishes (Order Siluriformes) are a particular diverse group of more than 3700 species within the superorder Ostariophysi (Fricke et al. 2021).Trichomycteridae is the second most diverse family of the order, comprising 349 valid species in nine subfamilies of which Trichomycterinae, shown to be monophyletic firstly by a morphological phylogenetic study (Datovo & Bockmann 2010) and more recently by molecular phylogenetic studies (Ochoa et al. 2017, 2020, Fernandez et al. 2021), is the most diverse with 253 valid species (Fricke et al. 2021).Trichomycterus Valenciennes, 1832 is the major genus with more than 160 species recognized (Froese & Pauly 2021), whereas the Andean genera Hatcheria Eigenmann, 1909 andBullockia Arratia, Chang, Menu-Marque &Rojas, 1978 are monotypic.The taxonomic relationships between Trichomycterus and Hatcheria have long been discussed (Arratia & Menu-Marque 1981, Fernandez et al. 2021).The genus Trichomycterus was erected in 1832 by Valenciennes (Cuvier & Valenciennes 1832), who in a following work described T. areolatus from the Maipo River basin, west of the Andes (Cuvier & Valenciennes 1846).In 1855, Girard described T. macraei collected near Uspallata, east of the Andes (Girard 1855), and then in 1909 Eigenmann created the genus Hatcheria and placed the specific epithet macraei within it (Eigenmann 1909).There were two other Hatcheria species described from Chilean waters, H. maldonadoi Eigenmann, 1927and H. bullocki Fowler, 1940, before Tchernavin (1944) put Hatcheria in synonymy with Trichomycterus, followed by De Buen who proposed it as a subgenus of Trichomycterus (De Buen 1958).Later, Ringuelet et al. (1967) mentioned five species of Hatcheria in Argentina; and Arratia et al. (1978) synonymized H. bullocki to H. maldonadoi including it in a new, monotypic genus for Chilean waters, Bullockia.Finally, Arratia & Menu-Marque (1981) synonymized all Argentinean Hatcheria species to H. macraei, leaving it as a monotypic genus.
The above taxonomic and nomenclatural arrangements have affected the identification of the geographic range for each of these species (Fig. 1a).The currently recognized distribution for T. areolatus well as their presence at one or both sides of the Andes have varied according to different authors.For T. areolatus, its southernmost distribution has been described at Puerto Varas (Pardo 2002), Abtao (Manriquez et al. 1988) and Chiloé Island at the X Región de los Lagos (Arratia 1981, Dyer 2000), as well as further south at General Carrera/Buenos Aires Lake in the XI Región de Aysén (Arratia et al. 1983).Moreover, it has also been identified in the eastern slope in the Cuyo region (Arratia & Menu-Marque 1981, Arratia et al. 1983).For H. macraei, its southernmost cis-Andean distribution is described at Blanco River (Unmack et al. 2012), as well as in the trans-Andean Baker and Aysén River basins (Arratia & Menu-Marque 1981, Campos et al. 1998, Unmack et al. 2012).Furthermore, recent records extended its trans-Andean northern distribution to the X Región de Los Lagos, where it was found in sympatry with T. areolatus in the Valdivia and Bueno River basins (Unmack et al. 2009a(Unmack et al. , 2012)).
Such taxonomic arrangements and variation in described range distributions seem to be due, at least in part, to the high morphological variability exhibited within these species, as well as similarities among them.Intraspecific differences in H. macraei were noted by Ringuelet et al. (1967) in his taxonomic key allowing the recognition of distinct populations, which were found to differ in body shape at different latitudes (Chiarello-Sosa et al. 2018).Likewise, Arratia et al. (1978) found meristic and morphological characters, and particularly those considering body relations, highly variable in B. maldonadoi.For T. areolatus, differences have been found between its northern (Choapa River) and central (Bío-bío River) populations 500 km apart (Pardo 2002), as well as between the latter and a southern population (Bueno River) separated by a similar distance (Colihueque et al. 2017).
Among the morphological characters used to distinguish between H. macraei and T. areolatus are caudal peduncle depth, dorsalray counts, position of anus, dorsal-fin length, and dorsal-fin shape (Arratia et al. 1978, Unmack et al. 2009b).However, these characters do not completely discriminate the species.For instance, Unmack et al. (2009b) found significant variation in caudal peduncle depth and position of anus in H. macraei.These authors used dorsalfin ray count as a nearly diagnostic character, with overall shape of the dorsal-fin as the only external character to discriminate both species.
The cause of such morphological variability is likely due to the aquatic environments (Senay et al. 2015).All three species occur in the rithronic zone of streams and rivers in areas with loose pebbles, gravel or sandy bottoms, substrates that allow individuals to bury themselves and avoid predators (Arratia 1983).Hatcheria macraei and T. areolatus prefer dark substrates, whereas B. maldonadoi prefer clearer bottoms (Arratia 1983).Bullockia maldonadoi is present in a few rivers with more restricted characteristics that are mainly pluvial and do not have tributaries (Arratia et al. 1978).In contrast, H. macraei and T. areolatus have broad environmental tolerances with a widespread presence in small headwater streams to low elevation rivers on a variety of substrates, especially gravel and rocks (Arratia et al. 1983, Unmack et al. 2009a, 2012).Moreover, it appears that habitat preference changes with age and can be somewhat flexible, allowing individuals to respond to seasonal river flow fluctuations.This may in part explain their widespread distribution (Arratia 1983).Adults of both species inhabit the benthic part of the rhitronal region of rivers and streams (Arratia 1983) and have been described as negatively phototactic and orienting themselves against the current (Ringuelet et al. 1967, Ringuelet 1975, Arratia 1976, Arratia & Menu-Marque 1981, Habit et al. 2005, Barriga et al. 2013).An important size-related habitat shift has been reported for H. macraei with positive allometric larval growth affecting mainly the head region, locomotion structures, the trunk region, and body robustness.Changes from positive allometric to isometric growth in H. macraei reflects the larva-juvenile transition (Barriga & Battini 2009).Juvenile and adult isometric growth has been described in T. areolatus (Habit et al. 2005).
Andean streams and rivers form many isolated basins that can promote population differences in low vagility organisms such as fishes (Hillman et al. 2014, Hancock & Hedrick 2018).Intrinsic riverine differences may prevail depending on the river basin (Pardo et al. 2005), which may be accentuated in large altitudinal and latitudinal ranges such as the geographical distribution of Andean trichomycterids.For instance, the northern distribution range of T. areolatus is characterized by short streams and rivers located about 4000 m above sea level (asl) with low flow and relatively high gradients, whereas southern rivers are relatively longer and at lower elevation, with relatively higher flows (Pardo et al. 2005).Such differences in riverine characteristics, including water velocity, have been shown to be associated with body shape differences in H. macraei (Chiarello-Sosa et al. 2018).
Even though morphological differentiation may or may not be present in closely related species, genetic divergence is expected as a result of genetic drift among species inhabiting isolated, environmentally different basins (Keeley et al. 2007).For instance, cryptic species sensu Mayr (i.e.species with poor morphological differentiation but that represent separate evolutionary lineages) are frequently observed in closely related taxa with broad distributions (Melo et al. 2016).In addition to phylogenetic inference, an analysis of geographic haplotype distributions may provide hints on the population structure and diversification between taxa.Moreover, species that have diverged very recently may prove difficult to distinguish by nuclear genes due to their relatively slow mutation rate (Hare 2001).On the contrary, the higher mutation rate of mitochondrial DNA makes it better fitted to resolve species-level and genus-level phylogenetic relationships among lineages (Rubinoff & Holland 2005).Being haploid and maternally inherited, mtDNA has one-quarter of the effective population size of nuclear genes, thus a mitochondrial haplotype tree can better resolve short internodes (Avise et al. 1987, Moore 1995) than other markers.
Here, we performed an integrative analysis to gain a more comprehensive knowledge of the limits of distribution of H. macraei, T. areolatus and B. maldonadoi.The goals of our work are: 1) to assess their relationship as well as their genetic variation and haplotype distribution in relation to hydrographic basins based on COI haplotypes -amplified for H. macraei in this study-and publicly available Cytb haplotypes for the three taxa generated by Unmack et al. (2009aUnmack et al. ( , 2012)); and 2) to analyze variation of body shape in H. macraei throughout its entire cis-Andean distribution in relation to Cytb genetic clades and river basins.

Collection of H. macraei specimens
We collected a total of 480 H. macraei individuals using seine net or electrofishing (Smith-Root backpack 24V, 600-900 V, 60 Hz, 6 ms standard pulse) in littorals of streams, rivers and lakes from 25 Argentinean localities in San Juan, Mendoza, Neuquén, Río Negro, Chubut, and Santa Cruz provinces, covering the whole cis-Andean distribution area (Fig. 1b, Table I).We sacrificed fish with an overdose of anesthesia (benzocaine 1:10000), and took digital images of the left side of each fish immediately to minimize parallax error for morphometric analysis.We preserved the samples in 96% ethanol to avoid the toxic formaldehyde, as well as to better preserve the DNA for molecular analysis.We sampled muscle tissue from 65 specimens (Supplementary Material -Table SI) collected from three localities at the Colorado River basin and one locality each at Negro and Yelcho River basins (Fig. 1b, Table I), before depositing them in Museo de La Plata (MLP), Facultad de Ciencias Naturales y Museo, Universidad Nacional de La Plata.We deposited the remainder specimens in the Department of Zoology at Centro Regional Universitario Bariloche, Universidad Nacional del Comahue, Argentina.We provide specimens deposit information in Table SI.Fish were caught and euthanized according to institutional guidelines for animal welfare and regulations detailed in Argentinean National Law No. 14346.Samplings were performed with the allowance of Parques

Phylogenetic assessment and genetic variability in H. macraei, T. areolatus and B. maldonadoi
We first tested the efectiveness of a nuclear gene to discriminate distinct evolutionary lineages on the three closely related trichomycterids here studied by running a Bayesian analysis for the RAG-1 nuclear gene utilizing 200 sequences available from GenBank (acc.n.JN186409-608).
Thereafter, we reassessed their relationships based on COI and Cytb mtDNA haplotypes.We extracted DNA from the 65 H. macraei muscle samples using the AccuPrep Genomic DNA Extraction kit (Bioneer, South Korea) and successfully amplified a 652-bp fragment of COI mtDNA gene from 55 of them (Table SI) by polymerase chain reaction (PCR) at the International Barcode of Life (iBOL) Argentinean reference Barcode Laboratory at the Museo Argentino de Ciencias Naturales (MACN) in Buenos Aires, Argentina.Amplifications were performed in a total volume of 12.5 μl consisting of 6.25 μl of 10% trehalose, 1.25 μl of 10X PCR buffer, 0.625 μl MgCl 2 (50 mM) 0.0625 μl of each dNTP (10 mM), 0.0625 μl of Taq DNA Polymerase (New England Biolabs), 2 μl of molecular grade water, 2 μl of DNA template, and 0.125 μl of each primer from the primer cocktail C_VF1LFt1-C_ VR1LRt1 (Ivanova et al. 2007).Cycling conditions consisted of an initial denaturation step of 94 °C for 2 min, followed by 35 cycles of 94 °C for 30 s, 52 °C for 40 s and 72 °C for 1 min, and a final extension at 72 °C for 10 min.We checked amplicons on 1.2% agarose gels before sending to sequencing on an ABI 3730XL capillary sequencer at the Canadian Centre for DNA Barcoding (CCDB) in Guelph, ON, Canada.We manually edited sequences with BioEdit v. 7.0.5.2 (Hall 1999) and made sequences available in the Barcode of Life Data Systems website (http:// www.boldsystems.org/).We identified nonredundant COI haplotypes using DnaSP v. 5.10 (Librado & Rozas 2009).We used MrBayes 3.1.2(Ronquist & Huelsenbeck 2003) to construct by Bayesian inference the phylogenetic relationship of the H. macraei haplotypes together with sequences from T. areolatus and B. maldonadoi (GenBank acc.n.KY857926, KY857963, and KY857964) previously published by Ochoa et al. (2017).We searched for COI sequences of all other southern Andean species that clustered together with H. macraei, T. areolatus and B. maldonadoi in clade E of Fernandez et al. 2021, however we found no sequences available.We then used sequences of Ituglanis parkoi (Miranda Ribeiro, 1944) (GenBank acc.n.KY857937) which clustered in sister clade D, and of T. brasiliensis Lütken, 1874 (GenBank acc.n.KY857993) that clustered in a more distant clade M from these authors's work, to root the tree.We identified the best-fit substitution model to be TIM1+I, nst = 6 and rates = equal, using the Akaike Information Criterion (AIC) with jModelTest 2.1.10(Posada 2008).We executed two independent runs totaling four chains, for 50 million MCMC generations, sampling every 1000 trees.We discarded the initial 25% of the resulting trees as burn-in, and constructed a 50% majority-rule consensus tree.We calculated pairwise genetic divergences among H. macraei haplotypes and T. areolatus, B. maldonadoi, I. parkoi and T. brasiliensis with PAUP* 4.0a (build 166) (Swofford 2003) using the above parameter estimates on the CIPRES Science Gateway (Miller et al. 2010).Finally, we constructed a haplotype network using the median joining algorithm developed by Bandelt et al. (1999), implemented in PopART v. 1.7 (Leigh & Bryant 2015), and edited it using the vector graphics editor Inkscape (https://www.inkscape.org).
For Cytb, we utilized the following haplotypes available from GenBank: H. macraei, n=73, acc.n.  2021) but found none available.Subsequently, we also included as outgroups sequences of Ituglanis parkoi (GenBank acc.n.KY858018/026) and T. brasiliensis (GenBank acc.n.KY858062).We identified the best-fit substitution model to be TIM3+I+G, nst = 6 and rates = gamma.We then ran a Bayesian inference analysis to calculate pairwise genetic divergences and constructed a haplotype network as described above using all the Cytb haplotypes.

Morphological variability in H. macraei
We carried out the morphological analysis on digital images from 447 individuals captured at 24 localities of Argentina (Fig. 1b, Table SI).We collected a total of 20 landmarks, four links, and five sliders (Fig. 2) on individuals with the software TpsDig version 2.05 (F.J. Rohlf, State University of New York, Stony Brook, NY).We aligned, rotated, translated and scaled landmarks through a Generalized Procrustes Analysis (GPA) using a consensus configuration (Rohlf & Slice 1990, Rohlf & Marcus 1993).We calculated partial and relative warps using Tpsrelw version 1.35, and visualized deformation grids using Relative Warps Analysis (RWA), thus using the consensus configuration as an average body shape (Cadrin 2000).We performed Discriminant Analysis (DA, Norusis 1986) employing partial warps and uniform coordinates (Rohlf 1996) in order to find body shape differences, considering capture sites, river basins that contain them, and Cytb clades obtained by Bayesian inference.We conducted further DA within each genetic clade considering capture sites as the grouping  (2023) 95(1) e20211007 8 | 24 variable.Basins were inappropriate as such due to the coexistence of two clades in some basins.

Phylogenetic assessment
Not surprisingly, we found the RAG-1 nuclear gene to be uninformative due to a lack of resolution (data not shown).On the contrary, mitochondrial genes were much more resolutive.We identified eight non-redundant haplotypes among the 55 COI sequences of H. macraei showing 19 variable sites (Table II).The COI Bayesian tree (Fig. 3) showed a well-supported subclade situating T. areolatus as the sister species of B. maldonadoi (posterior probability [PP] = 0.92).With the exception of haplotype 7 (Rosario Lake, Yelcho River basin, a Pacific drainage), relationships among the H. macraei COI haplotypes were shallow, although they did reveal a geographic pattern clustering southernmost haplotypes 6, 7, and 8 along T. areolatus and B. maldonadoi (PP = 0.94).The Colorado River basin haplotypes 1-5 remained basal to that subclade (Fig. 3).
The Cytb Bayesian tree resulted in four well-supported major clades, all in a polytomy (Fig. 4).The "Bullockia spp.Clade" presented three subclades.The first one clustered the northernmost haplotypes from the Limarí, Choapa, Petorca and Aconcagua River basins (previously assigned to T. areolatus clade "A" by Unmack et al. 2009a).The second subclade "B.maldonadoi", clustered all B. maldonadoi haplotypes including those from its type locality (i.e.Andalién River; Arratia et al. 1978).The third subclade clustered three haplotypes from the Reloca River along two haplotypes from the Maipo River (previously assigned to T. areolatus clade "B" by Unmack et al. 2009a).In light of the pairwise genetic distances (see below), we labelled the subclades closely-related to the one of B. maldonadoi as "B.m. ssp 1" and "B.m. ssp 2" until further studies are conducted to resolve their taxonomic relationship to B. maldonadoi.The second, "T.areolatus Clade", clustered many of T. areolatus haplotypes from Imperial and Toltén River basins in central Chile, along most of those from northern freshwater systems up to the Maipo River.It was composed of three subclades corresponding to clades "C", "D", and "E" sensu Unmack et al. (2009a).The third clade contained a minor subset of T. areolatus haplotypes from three central rivers: Imperial, Toltén and Valdivia (clade "F" sensu Unmack et al. 2009a).We here refer to it as "T.areolatus ssp.
The haplotype network for the Cytb mtDNA region (Fig. 5) exhibited high resolution as expected, as the haplotypes employed were representative of the whole distribution range of the southern Andean trichomycterids here considered.This network showed a deep genetic structure with a likely ancestral haplogroup common to all present haplogroups.In particular for H. macraei, the Cytb haplotype network showed three haplogroups coincident with the Bayesian tree for this marker: Atlantic, Pacific and Atlantic-Pacific haplogroups.On the contrary, the COI haplotype network for H. macraei presented a poor resolution as the sampling was limited to a few populations in their cis-Andean range.There must be undetected haplotypes that if included would improve the outcome of this network.

Geometric morphometrics
The first two RWs of the GMA explained 31.9% of the variance (RW1= 22.5% and RW2= 9.4%), indicating a mostly isodiametric shape of the morphological data cloud in the hyperspace.The body shape variation explained by RW1 ranges between individuals with a bigger head, shorter trunk, a shorter basis of the dorsal-fin and a higher caudal peduncle in the positive semi-axis, and individuals with a smaller head, longer trunk, a longer base of the dorsal-fin and a lower caudal peduncle, in the negative semiaxis.The body shape variation explained by RW2, ranges between individuals with a bigger head, a shorter base of the dorsal-fin, and higher caudal peduncle in the negative semi-axis, and individuals with a smaller head, a longer basis of the dorsal-fin and lower caudal peduncle in the positive semi-axis (Fig. 6).At odds with the RWs, the DA of the partial warps and uniform coordinates, grouped data successively by capture sites, rivers basins, and Cytb genetic clades for H. macraei, and showed the latter to give better discrimination with two significant discriminating functions explaining 100% of the variance, and 95.7% of original grouped cases correctly classified (Fig. 7, Table V).
We subsequently processed the original database separately within each of the three  grouped cases correctly classified.However, a significant cubic regression (P< 0.0001) of the latter two DF1 and DF2 with body size (total length) for the latter clade implies that these differences must be considered with caution.All files and photographs employed in the morphometric analysis are available at https:// ri.conicet.gov.ar/handle/11336/151708.

DISCUSSION
In this integrative study, we analyzed the morphologic variation of H. macraei and the genetic variation of H. macraei, T. areolatus and B. maldonadoi, in relation to hydrographic basins to help further our knowledge of their limits of distribution.Our results unveiled new geographic boundaries between B. maldonadoi and the northernmost T. areolatus populations, as well as between H. macraei and the southernmost T. areolatus populations.Recent advances in knowledge of catfishes in South America (Katz et al. 2018, Ochoa et al. 2017, 2020, Costa et al. 2020, Fernandez et al. 2021) allowed us to examine the relationships among lineages of these three southern Andean trichomycterids.By using all Cytb haplotypes available for them in a single phylogenetic inference, our results revealed some lineages previously assigned to T. areolatus actually have a closer relationship to either B. maldonadoi or H. macraei.Moreover, our results confirmed the monophyly of Hatcheria, which clustered all of its lineages together, and of Bullockia, by which it was also confirmed the monophyly of B. maldonadoi and its distinctiveness by presenting the relatively longest branch (Fig. 4), coincident with results in Fernandez et al. (2021).Our analyses resolved Bullockia as the sister group of T. areolatus, contrary to Fernandez et al. (2021) who found it to be H. macraei.First, our Bayesian inference tree for the COI marker showed B. maldonadoi and T. areolatus clustering together with high support (Fig. 3, PP=0.92).Second, the Cytb pairwise genetic distance of T. areolatus from Bullockia spp. was lower than that from H. macraei, being 5.2% and 6.1% respectively (Table IVa).Third, the Cytb haplotype network showed Bullockia haplogroups more closely related to Trichomycterus than to Hatcheria (Fig. 4).Finally, and contrary to the current knowledge, we found that the geographic distribution of populations belonging to Bullockia spp.and T. areolatus each span an extent of about 6° of the meridian, with about half of that range overlapped.
Our results amount to the emerging evidence for geographically circumscribed subclades of Trichomycterus species as described by Fernandez et al. (2021), and provide hints this may also occur in Bullockia, which its fragmented distribution is reflected by the high divergence among its lineages (4.1%;Table IVa).Across the vast range of Bullockia, from the Limarí River at 30°40'S south to the Ramadillas River at 37°18'S (Fig. 5), B. maldonadoi showed a disjunct distribution in the very south coincident with was previously described by Arratia et al. (1978), ranging from the Itata River at 37°09'S to the Ramadillas River (Fig. 5).Both B. maldonadoi ssp. 1 and B. maldonadoi ssp. 2 presented relatively high genetic distances of 4.2% from B. maldonadoi, with even a significant distance of 3.8% between themselves (Table IVb).This brings attention to the current status of Bullockia as a monotypic genus, as it can be hypothesized B. maldonadoi may have sister species north from the Itata River.We could argue that the northern Bullockia populations in Limarí, Choapa, Petorca and Aconcagua River basins belong to a cryptic species, or to a subspecies of B. maldonadoi, in the same sense we could argue they are Operational Taxonomic Units (OTUs).Any of these alternatives are just pragmatic definitions to group individuals by similarity (Sokal & Crovello 1970;Mayr 1982).Therefore, we highlight the necessity of further investigation to describe patterns of character distributions of Bullockia spp.subclades.
To date, the geographic distribution of T. areolatus was considered to be between 30°40'S (Limarí River) and 42º 22'S in Chiloé Island (Pardo 2002, Unmack et al. 2009a, b).However, our results showed it to range about half of that, between 33°S in Maipo River basin -which contains its type locality (Girard 1855) -and 39°S at the Toltén River.This river is considered the northern limit of a geological province, the Patagonian Cordillera, characterized by a batholitic belt that extends southward to Cape Horn Islands (Ramos & Gighlione 2008).Besides its peculiar geological significance, the Toltén River also holds a strong climatological meaning as it is located at the estimated maximum extent of the major ice sheet during the last glacial maximum (LGM) (Rabassa et al. 2011).Moreover, our results showing Toltén River as the actual southern boundary of T. areolatus distribution agrees with the geographic limit Dyer (2000) established between the South-Central and the Southern ichthyogeographic subprovinces.
The great morphological variation among H. macraei populations was consistent with the intra-specific variation described in studies by Arratia & Menu-Marque (1981) and Chiarello-Sosa et al. (2018).This was reflected mainly in the head size, trunk length, caudal peduncle height, and dorsal-fin length.Different microhabitat river conditions (e.g.water velocity) can affect body shape, mostly body height and fin size, as shown by the work of Chiarello-Sosa et al. (2018).These authors found a relationship between stream flow and H. macraei body shape, in which the faster the water flow, the narrow the caudal peduncle is as well as the longer the base of the dorsal-fin.In this sense, the morphological variation we described in H. macraei was congruent with river basins and Cytb genetic clades (Figs. 7 and 8; Table V), resulting in three well defined morphological groups.One group was formed by the northern Colorado River basin within the Atlantic Clade, another one was extended southward from the Negro River basin within the Atlantic-Pacific Clade, and the other group was restricted to two populations separated by more than 5 latitudinal degrees within the Pacific Clade (Lake Cholila in the head waters of the Yelcho River; and Blanco River -its southernmost record, presently a closed basin).
Immediately southward from the Toltén River, the identified southern boundary of T. areolatus, we confirmed the Valdivia River as the northern boundary of the trans-Andean distribution of H. macraei as previously reported by Unmack et al. (2009aUnmack et al. ( , 2012)).Moreover, our data showed that the remaining trichomycterid populations analysed to the south in Bueno, Llico and Maullín Rivers, as well as those from Butalcura River in Chiloé Island, which were previously assigned to T. areolatus in Unmack et al. (2009a), are actually best aligned with H. macraei, thus expanding the southern boundary of the trans-Andean distribution of this species (Fig. 5).Furthermore, the position of haplotypes from these four populations at the base of the Atlantic-Pacific subclade of H. macraei (Fig. 4), suggests that these rivers may have been refugia during the later glacial periods.Evidence of glacial refugia has been found in Coastal Chile and southern flanks of the Andes (Premoli et al. 2000), and populations located outside ice limits seem to have been isolated during the glacial period (Premoli et al. 2002).Following recent post-glacial periods, the subsequent deglaciation and drainage reversals (Unmack et al. 2012) may have allowed H. macraei to colonize freshwater systems at the eastern slope of the Andes.This scenario can also be interpreted from the shallow nodes that characterize the Negro, Chubut and Deseado River basins in the Atlantic-Pacific subclade (Fig. 4), as well as their star-like shape in the Cytb haplotype network, showing them as recent expansions (Fig. 5).
The southern Andean trichomycterids here studied diversified in drainage networks characterized by strong flow rates that formed during the Andes uplift, which generated the vast extension of gravel mantels that cover most of Patagonia (Martínez & Kutschker 2011).Lineages within each species have likely evolved in allopatry as a consequence of the history of glacial cycles in southern South America during the Cenozoic (Rabassa et al. 2011).Therefore, the finding of small lineages such as the T. areolatus ssp.over an LGM boundary is not surprising since this type of area, which likely harbored refugia during glacial cycles, may indeed concentrate many endemisms (Dyer 2000).This lineage presented a genetic distance of 2.8% from the major T. areolatus Clade, doubling the within-clade distance of the latter (Table IVa).Taking into consideration that T. areolatus ssp.also presented a low within-group divergence of 0.3% (Table IVa), the fact that they live in sympatry with T. areolatus in Imperial and Toltén River basins, and with H. macraei in the Valdivia River, it would be prudent to consider this small lineage as a subspecies of T. areolatus until further evidence is collected.We hypothesize this divergent lineage may have originated as a consequence of isolation in a glacial refugium in the Valdivia River.The star-like shape of the haplotype network for this species, with a high number of haplotypes, and its most common haplotype the one of a few present in the Imperial and Toltén River basins in the north (Fig. 5), seem to support this.Fossil evidence suggests that Siluriforms were widely distributed in southern South America until at least the late Miocene, with a warmer climate than today (Cione et al. 2005, Dozo et al. 2010, Azpelicueta & Cione 2016).Starting in early Pliocene, a regional uplift of the Andes promoted the eastward aridization of Patagonia (Ramos & Ghiglione 2008).Although warm climate could have extended until late Pliocene (Cione et al. 2005), glaciation cycles (Rabassa et al. 2011) likely fragmented further their distribution.Several putative refugia for plant and vertebrates seem to have existed in southern South America during the climate oscillations that characterized the Pleistocene (Sérsic et al. 2011), which may have isolated to and reunited populations several times.Therefore, the Andean uplift followed by glacial cycles may have promoted the origin of many endemisms in trichomycterids, which may not necessarily involve morphological change.Morphologically static cladogenetic events (Bickford et al. 2007) can result in morphologically indistinguishable species which are taxonomically classified as being one, i.e. cryptic species (Mayr 1942).There are many instances in recent years of cryptic species described in Neotropical fishes (e.g.Martin & Bermingham 2000, García-Dávila et al. 2013, Gomes et al. 2015, Melo et al. 2016, Guimarães et al. 2020, Pereira et al. 2021), as well as in other regions of the world (e.g.Africa: Jirsová et al. 2019, Asia: Shao et al. 2021).Looking at the population structure of marinedescendant species, that settled well after the presence of the Siluriformes in Patagonia, can allow us to compare groups of taxa at the present time.The coincidence observed here between morphological grouping and genetic clades at the level of mitochondrial markers in H. macraei, also demonstrated in populations of T. areolatus (Pardo et al. 2005), does not occur in marine related Patagonian fishes.Both Percichthys trucha (Valenciennes, 1833) and Odontesthes hatcheri (Eigenmann, 1909) show a major degree of morphological variation (Crichigno et al. 2012, 2013, 2016, Conte-Grand et al. 2015).However, mitochondrial markers failed to reveal genetic structuring between populations for these species (Ruzzante et al. 2006, 2011, Conte-Grand et al. 2015, Rueda et al. 2017).Therefore, we suggest the older cladogenetic events that shaped freshwater assemblage in the southern Andean trichomycterids (Barber et al. 2011) prevailed over the events that homogenized the populations of P. trucha and O. hatcheri (Ruzzante et al. 2006, 2011, Crichigno et al. 2014, Conte-Grand et al. 2015, Rueda et al. 2017).Not surprisingly, another Patagonian species that showed an appreciable genetic structure at mitochondrial level, and evidenced refugia into the glaciated area, is Galaxias platei Steindachner, 1898(Zemlak et al. 2008).Its presence in South America has been dated not just before the Andes upraise in the Miocene (Ruzzante et al. 2008;Zemlak et al. 2008), but even before the separation of Australia and South America, in the Oligocene (Burridge et al. 2012).In light of our results, the historical taxonomic and nomenclatural arrangements of the relict southern Andean trichomycterids here investigated seem to have originated by their great intra-specific morphological variation, while at the same time by their high inter-specific morphological similarity.Both can be interpreted as consequences of sharing a recent common ancestor, as well as the great variety of microhabitat where these species live in lotic environments, from mountain rocky streams with rapid waters to slow-habitats with sandy bottoms (Chiarello-Sosa et al. 2018).All in all, we confirmed both the great morphologic variation of H. macraei and the high genetic variation of H. macraei, T. areolatus and B. maldonadoi, to be congruent with hydrographic basins.We also revealed new boundaries to the currently known trans-Andean distribution of these relict trichomycterids and highlight the need of further integrative studies to continue improving our knowledge of the southern Andean trichomycterid diversity.

Figure 3 .
Figure 3. Bayesian phylogenetic tree based on eight Hatcheria macraei haplotypes obtained by amplification of the COI mtDNA region, as well as two haplotypes from Trichomycterus areolatus and one from Bullockia maldonadoi available at GenBank.Haplotypes from T. brasiliensis and Ituglanis parkoi were used as outgroups.GenBank acc.n. of sequences for B. maldonadoi, T. areolatus, and outgroup taxa are indicated in their respective terminals.Bayesian posterior probabilities are shown on nodes.

Figure 4 .
Figure 4. Bayesian phylogenetic tree based on 73, 125, and 12 publicly available haplotypes of the Cytb mtDNA region for Hatcheria macraei, Trichomycterus areolatus and Bullockia maldonadoi, respectively.Haplotypes from T. brasiliensis and Ituglanis parkoi were used as outgroups.Bayesian posterior probabilities are shown on major nodes.Terminals are represented by haplotypes with their corresponding GenBank accession numbers and the river basins where are present (see Fig. 5 for color coding).Clades and subclades are indicated to the right of the tree.Letters in brackets A-G correspond to the seven major clades described for T. areolatus by Unmack et al. (2009a).Abbreviations in brackets (BC: Big Clade; B-Ch: Baker-Cholila; Co: Colorado) correspond to Cytb clades described for H. macraei by Unmack et al. (2012).

Figure 5 .
Figure 5. Median joining networks constructed with eight haplotypes of the COI mtDNA region found in Hatcheria macraei individuals in the present study (upper right); and with 73, 125, and 12 publicly available Cytb haplotypes from H. macraei, Trichomycterus areolatus and Bullockia maldonadoi, respectively (left side).Nodes represent unique haplotypes scaled according to their frequency, with colors matching with the basin where they are present (for location of rivers see also Fig. 1a).Black nodes represent inferred unsampled or ancestral haplotypes.Branch lengths are proportional to the number of nucleotide mutations between nodes.Single nucleotide mutations between nodes are indicated with bars crossing the branches if mutations <= 5, or with numbers in boxes if mutations > 5. Letters in brackets A-G correspond to the seven major Cytb clades described for T. areolatus by Unmack et al. (2009a).Abbreviations in brackets (BC: Big Clade; B-Ch: Baker-Cholila; Co: Colorado) correspond to Cytb clades described for H. macraei by Unmack et al. (2012).

Table I .
Sampling localities for Hatcheria macraei with abbreviation codes employed in this study, their geographic coordinates and drainage basins.Number of individuals captured at each sampling site (N S ), number of samples successfully amplified for COI (N COI ) and number of specimens included in the morphometric analysis (N M ).

Table II .
Haplotypes and nucleotide position of variable sites found in COI sequences amplified from 55 individuals of Hatcheria macraei.
Cytb clades found Bullockia spp. as the most divergent, with

Table III .
Pairwise genetic distances corrected under the best-fit GTR substitution model calculated for COI, comparing the eight Hatcheria macraei haplotypes described in this study and haplotypes available from GenBank for Trichomycterus areolatus and Bullockia maldonadoi.Comparisons with the more distantly related species Ituglanis parkoi and T. brasiliensis are included.

Table IV .
Pairwise genetic distances corrected under the best-fit GTR substitution model calculated for Cytb comparing: a) major Cytb clades; b) subclades within the Bullockia spp.clade; c) subclades within the Trichomycterus areolatus clade; and d) subclades within the Hatcheria macraei clade.Mean genetic distances are indicated in bold, with corresponding ranges in parentheses below.

Table V .
Discriminant analyses of the partial warps and uniform coordinates obtained from geometric morphometrics analysis for Hatcheria macraei (N= 447), by capture sites, river basins, and Cytb genetic clades.