Drainage Network Morphology Influences Population Structure and Gene Flow of the Andean Water Frog Telmatobius pefauri (Anura: Telmatobiidae) of the Atacama Desert, Northern Chile

© 2023 Academia Sinica, Taiwan Open Access Drainage Network Morphology Influences Population Structure and Gene Flow of the Andean Water Frog Telmatobius pefauri (Anura: Telmatobiidae) of the Atacama Desert, Northern Chile Pablo Fibla1, Paola A. Sáez1,2 , Franco Cruz-Jofré1,3 , and Marco A. Méndez1,2,4,* 1Laboratorio de Genética y Evolución, Facultad de Ciencias, Universidad de Chile. *Correspondence: E-mail: [email protected] (Méndez) E-mail: [email protected] (Fibla); [email protected] (Sáez); [email protected] (Cruz-Jofré) 2Center of Applied Ecology and Sustainability (CAPES) 3Facultad de Recursos Naturales y Medicina Veterinaria, Universidad Santo Tomás, Chile 4Instituto de Ecología y Biodiversidad (IEB), Santiago, Chile Received 6 December 2022 / Accepted 26 June 2023 / Published 1 September 2023 Communicated by Teng-Chiu Lin Desert aquatic species tend to show isolated and disconnected populations due to the fragmented nature of their environment; however, the morphology of the hydrographic basins, added to humid climatic conditions, can allow dispersion between populations in a desert environment. The aim of this study was to examine the influence of drainage morphology on the phylogeographic structure and gene flow (using a fragment of the mitochondrial control region and seven microsatellite markers) of an endemic taxon of the Andean Precordillera in the Atacama Desert, the aquatic frog species Telmatobius pefauri. We detected three genetic clusters, one cluster present in the Lluta basin and two clusters in the Azapa basin. The results suggest that the genetic structure of T. pefauri is influenced by the morphology of the drainage network formed by the Lluta and Azapa basins: localities present in the same drainage, Tignamar River, were less differentiated and showed higher gene flow levels among them than to their conspecifics belonging to the other drainage in the same basin, Seco River, and those belonging to the other basin, Lluta basin. Gene flow patterns and genetic structure to populations Atacama Andean aquatic taxa would be influenced by basin morphology, with dispersion being stimulated in dendritic hydrological systems, and eventually by humid climatic (regional) events. Key words: Water frog, Andean Precordillera, Atacama Desert, Dendritic hydrographic network, Stream Hierarchy Model Citation : Fibla P, Sáez PA, Cruz-Jofré F, Méndez MA. 2023. Drainage network morphology influences population structure and gene flow of the Andean water frog Telmatobius pefauri (Anura: Telmatobiidae) of the Atacama Desert, Northern Chile. Zool Stud 62:44. doi:10.6620/ZS.2023.62-44. BACKGROUND For aquatic species, desert environments greatly restrict the gene flow between disconnected bodies of water, generating marked phylogeographic structure and high levels of differentiation between populations, which is exacerbated in species that have a low dispersal capacity (Meffe and Vrijenhoek 1988; Seidel et al. 2009; Morales et al. 2011; Guzik et al. 2012; Murphy et al. 2012 2013). However, there are cases in which drainage morphology and humid climatic episodes can promote gene flow between populations of aquatic taxa present in desert systems (Meffe and Vrijenhoek 1988; Hughes et al. 2009; Collado and Méndez 2013; Vila et al. 2013; Cruz-Jofré et al. 2016). The isolated nature of creeks or desert springs, Zoological Studies 62:44 (2023) doi:10.6620/ZS.2023.62-44 1 © 2023 Academia Sinica, Taiwan combined with the low dispersal capacity of aquatic organisms in an arid landscape, can result in complex gene flow patterns, where the degree of connectivity between creeks would explain to a greater extent a given phylogeographic structure (Unmack et al. 2013). Meffe and Vrijenhoek (1988) proposed two models that can describe the connectivity of taxa that inhabit desert wetlands and make predictions about their diversity and genetic differentiation between populations: the Death Valley Model and the Stream Hierarchy Model. The Death Valley Model (hereafter DVM) refers to remnant populations without terrestrial dispersion that are completely restricted to their native, isolated springs. Strong isolation in combination with small local population size is expected to cause local genetic drift to dominate population genetic structure in the absence of the homogenizing force of dispersal (Hughes et al. 2009). Under the DVM, populations should show low genetic diversities and high levels of differentiation. On the other hand, the Stream Hierarchy Model (hereafter SHM) is representative of a dendritic system, with gene flow between populations facilitated by direct hydrological connectivity. Under this model, the degree of isolation should follow a hierarchical pattern of genetic structure that reflects the nested arrangement of subcatchments within catchments within major basins of the stream network; genetic variance should partition significantly among drainage basins at any spatial scale at which basins can be defined within the stream network (Hughes et al. 2013). The Andean zone of northern Chile (dry Central Andes) is embedded in the South American Arid Diagonal (Villagrán and Hinojosa 1997) and forms an important part of the Atacama Desert. The main geological divisions of this zone are (from west to east): the Intermediate Depression (“Atacama”), the Precordillera and the western Andean flank, and the Altiplano (Placzek et al. 2009), which reaches an average altitude of 4,000 masl (Garzione et al. 2008). During the austral summer (DJF), in the Atacama Desert, humid episodes occur along the Western Andean Cordillera and the Altiplano as a result of the South American Summer Monsoon (SASM), when strong upper-level easterly winds enhance moisture transport from Amazonia, creating saturation during uplift within deep convection cells (Garreaud et al. 2003). As a consequence of the source of humidity from the east, a “rain shadow” effect is generated over the Western Cordillera and the Atacama Desert, where the average annual precipitation decreases rapidly in the altitudinal gradient, with more than 300 mm per year at 5,000 masl to less than 20 mm per year at 2,300 masl. Below 2,300 masl, associated with the Central Valley (or Intermediate Depression), is an area of extreme hyperaridity in which the average annual precipitation is less than 1 mm per year (Houston and Hartley 2003). The bodies of water in this area that are above 2,300 masl are often geographically isolated without hydrological connections and separated by long distances (Vila et al. 2006). From a hydrological point of view, this region is fragmented into a series of hydrographic basins that are partially or totally isolated from each other and delimited by irregular reliefs, mountain ranges and even volcanoes. The majority of the basins in northern Chile are endorheic, mainly associated with the high plateau zone and the Intermediate Depression, and there are only a few exorheic basins in the Precordillera in the extreme north of Chile (Mortimer 1980). The high degree of geographic isolation of the wetlands present in the Andean zone of northern Chile and the presence of geomorphological barriers between the hydrographic basins would have stimulated the origin of a high number of aquatic endemism in this region, especially in systematic groups such as amphibians and fishes, as well as in aquatic invertebrates (Vila et al. 2006 2013; Collado et al. 2011; Sáez et al. 2014). For several of the endemic taxa of northern Chile, an evolutionary scenario has been proposed in which the populations of different freshwater species were progressively fragmented along with the aquatic environment at the end of the Pleistocene and the beginning of the Holocene, which has led to the suggestion that population divergence/ differentiation events in this area would be mainly associated with the fragmentation of ancestral water bodies (Collado et al. 2011; Vila et al. 2013; Sáez et al. 2014). This hypothesis is mainly supported by aquatic organisms present in highly fragmented highland systems (such as in disconnected springs in the Ascotán salt pan) in which high levels of genetic differentiation, marked phylogeographic structure and restricted gene flow between units has been observed (Morales et al. 2011; Cruz-Jofré et al. 2016; Valladares et al. 2018). However, there are populations of aquatic organisms (amphibians and fish) in the southern Altiplano that inhabit an extensive drainage network and that show high hydrological connectivity (e.g., Orestias and Telmatobius populations in the Lauca River), in which low genetic divergence has been detected between populations (Victoriano et al. 2015; Guerrero-Jiménez et al. 2017). This suggests that the drainage pattern and morphology of the basin could have influenced the genetic structures of the species present in northern Chile, inviting the possibility of detecting gene flow (including founder events) between freshwater populations in systems that show hydrological connectivity to lower elevations (relative to the Altiplano) within the Atacama Desert. page 2 of 12Zoological Studies 62:44 (2023) © 2023 Academia Sinica, Taiwan High Andean environments are the habitat of one of the most diverse groups of amphibians in South America, the genus Telmatobius Wiegmann, 1834. The species of the genus Telmatobius present in Chile have strictly aquatic habits and inhabit freshwater systems that show different drainage patterns, either in the Precordillera or in the Altiplano, with some species associated with the Atacama Desert (Veloso and Trueb 1976; Veloso et al. 1982; Formas et al. 1999 2003 2005 2006). Due to their limited vagility and high philopatry, amphibians generally have highly genetically structured populations over short geographic distances, and thus often retain signals of historical events that have led to their current distributions (Vences and Wake 2007; Zeisset and Beebee 2008). The amphibian Telmatobius pefauri inhabits the pre-Andean semidesert zone of extreme northern Chile. Populations of this species are distributed between 3,000 and 3,500 masl in two of the main exorheic drainages of the region, the Lluta and Azapa basins (Fibla et al. 2017). The Lluta and Azapa Rivers form the main branches of a trellis network drainage pattern, the catchments of which are constituted by two dendritic hydrographic networks separated by the Cordon Quevilque mountain ridge (Fig. 1; García and Heraíl 2005). The known populations of T. pefauri are present in high altitude tributaries in the middle and higher zones of the Lluta and Azapa Basins, respectively. Most populations of T. pefauri are present Fig. 1. Study area, distribution of Telmatobius pefauri. Localities, 1: Socoroma (Socoroma River); 2: Murmuntani; 3: Copaquilla; 4: Chapiquiña; 5: Belén; 6: Lupica; 7: Saxamar. Localities 2 and 3 belong to the Seco River drainage; localities 4–7 belong to the Tignamar River drainage. Basin limits are indicated with dashed lines. The inset map shows the study area (highlighted by a red box) in relation to South America. SAAD = South American Arid Diagonal. page 3 of 12Zoological Studies 62:44 (2023) © 2023 Academia Sinica, Taiwan in the Azapa basin, in localities belonging to either the Seco River or the Tignamar River drainages; just one population of this species is known in the Lluta basin. These characteristics make T. pefauri an excellent model to evaluate the influence of drainage morphology and climatic events on the phylogeographic structure and gene flow of the endemic taxa of the Andean Precordillera and the Atacama Desert. Given the basin morphology and drainage network pattern of T. pefauri distribution, the climatic history of the study area and the expected null ability of this species to disperse out of water bodies, we expect to find that populations of T. pefauri meet the prediction of the SHM: the existence of high genetic differentiation between populations belonging to different drainages in the same network (i.e., Seco and Tignamar Rivers) and/or to different drainage networks (Lluta and Tignamar basins), and greater gene flow levels between populations belonging to the same drainage and/or drainage network. To study this, we applied a phylogeographic approach and conducted coalescence analysis based on a fragment of the mitochondrial control region and microsatellite markers. MATERIALS AND METHODS Specimen sampling and molecular analysis During January 2017 and March 2018, biological samples of T. pefauri specimens (larvae, juveniles and adults) were obtained from locations within the species’ known distribution (Fig. 1). Sampled localities correspond to: Socoroma, Murmuntani, Copaquilla, Chapiquiña, Belén, Lupica, Saxamar; coordinates are described in Fibla et al. (2017), except for Chapiquiña (18°23'27.40"S; 69°32'25.57"W). Frogs were captured under the vegetation and among the rocks of the creeks using fishing nets. After capture, specimens were anaesthetized using benzocaine ([ethyl 4-aminobenzoate]; 300 mg/L) diluted in water from the same source where they were captured. A small piece (< 0.1 ml) of interdigital membrane (adults and juveniles) or tail membrane (larvae) was extracted from each specimen. Tissue samples were stored in 2 ml tubes filled with 100% ethanol. Specimens were released at the same collection site immediately after full recovery from anaesthesia. Total DNA was isolated using the salt extraction method (modified from Jowett 1986). Verification of DNA extraction and quality was performed by agarose gel electrophoresis (2%) with GelRedTM staining (Biotium) visualized with UV light. DNA quantification was performed on a NanoDrop LiteTM spectrophotometry kit (Thermo Scientific). For genetic analysis, a fragment of the mitochondrial control region and microsatellite (SSR, simple sequence repeat) markers were amplified. The primers used for the amplification (PCR) of the control region fragment were Tchu_cytb (1103) (forward; CAACAATCGGAGCACTAGA) and Tchu_ dloop (1197) (reverse; CCTAGCTCCTGACTTCTT). The amplification conditions (PCR) were as follows: 3.0 mM MgCl2, 0.4 mM dNTPs, 0.15 mM primers, 1.0 U Taq and 60 ng of total DNA in a final volume of 30 μL. The thermal profile was 3 min at 94°C for initial denaturation, followed by 40 cycles of 30 s at 94°C, 45 s at 54°C–58°C for annealing, 80 s at 72°C, and lastly, 10 min at 72°C in the final extension phase. The sequences obtained (approx. 993 bp) were aligned and manually edited in the program BioEdit v.7.2.0 (Hall 1999). The DNA sequence matrix used in the analyses was constructed using the ClustalW algorithm for multiple alignments. Six individuals collected in Socoroma (Lluta basin) showed DNA sequences from another Telmatobius species, presumably Telmatobius marmoratus. Therefore, they were excluded from further analyses. In the case of microsatellite markers, eight markers of the roughly 30 primer pairs tested showed polymorphism, most of which were selected from Fabres et al. (2018). Selected markers were amplified following the amplification conditions (PCR) and labelling protocol detailed by these authors. The genotypes of the microsatellite loci were annotated manually using GeneMarker software (Softgenetics, State College, PA). All loci within each population were checked for the presence of genotyping errors, null alleles, large allele dropout and stuttering with microchecker v.2.2.3 (Van Oosterhout et al. 2004) using Bonferroni-adjusted 95% confidence intervals (DunnSidak) derived from 10,000 Monte Carlo simulations. Linkage disequilibrium (LD) between pairs of loci was calculated with GENEPOP v.4.4.2 (Rousset 2008) using the Markov chain method (10,000 de-memorization steps, 1,000 batches and 5,000 iterations/batch) and a sequential Bonferroni correction. Genetic haplotype diversity and intraspecific genealogies For mitochondrial sequences, the levels of genetic polymorphisms were estimated by determining the number of haplotypes (K), the number of polymorphic sites (S), haplotype diversity (H), the average number of differences between pairs of sequences (Π) and the nucleotide diversity (p) for each locality using DnaSP v.5.0 (Librado and Rozas 2009). A network of page 4 of 12Zoological Studies 62:44 (2023) © 2023 Academia Sinica, Taiwan haplotypes was reconstructed using the median-joining algorithm (Bandelt et al. 1999) implemented in PopART (http://popart.otago.ac.nz/index.shtml; Leigh and Bryant 2015). The diversity by locality of the SSR markers was evaluated in terms of the average number of alleles per locus (NA), observed heterozygosity (HO) and expected heterozygosity (HE) using GENETIX v. 4.05.2 (Belkhir 2004). Deviations from Hardy-Weinberg equilibrium were evaluated by calculating the inbreeding coefficient (FIS) by locality using 10,000 permutations to obtain their significance. Population genetic structure The degree of genetic differentiation between populations was determined by estimating the pairwise FST between localities in Arlequin v.3.5 (Excoffier and Lischer 2010) with mitochondrial data and in GENETIX with microsatellites, using 10,000 permutations to evaluate their significance. The Bonferroni correction was applied to adjust the p-values of the multiple comparisons. The GENELAND version 2.0 package (Guillot et al. 2005) was used to estimate the number of panmictic clusters and locate their spatial limits based on mitochondrial and microsatellite data. For this analysis, 15 independent MCMC runs were used, each with 5,000,000 iterations sampled every 500 steps and considering a burn-in of 25%, evaluating between one and seven clusters (CTmin = 1 and CTmax = 7). To corroborate Geneland results without considering spatial information (since hydrologic distance could be different than geographic distance among T. pefauri’s populations), a Principal Components Analysis was performed using the Adegenet package (Jombart 2008) in R environment (R Core Team 2012). Gene flow and historical population sizes To estimate the levels of historical gene flow between the inferred genetic clusters and their population sizes, a Coalescent Bayesian approximation was implemented in MIGRATE-N 4.2.1 (Beerli and Palczewski 2010). Migrate analysis was performed using Markov Chain Monte Carlo (MCMC) and consisted of a single long chain with 10 replicates, 1,000,000 values of parameters visited (genealogies) and a total of 10,000 recorded steps (saved), considering a burn-in of 1,000 trees per chain. Heating was performed under a static scheme with four temperatures (1, 1.5, 3 and 100,000). Ten independent runs were performed under the HKY substitution model (suggested by JmodelTest 2.0; Darriba et al. 2012). The Theta adjusted population sizes (θ = NE*μ, where NE is the historical effective size and μ is the mutation rate per generation) and the adjusted historical migration rate M (M = m/μ, where m is the immigration rate per generation) were estimated from the best run, which was selected in terms of its Bezier Approximated log-marginal likelihood score. The number of effective migrants per generation was calculated as NE*m = θ*M. To evaluate the recent gene flow (last three generations) between localities, BayesAss 3.0 (Wilson and Rannala 2003) was used using SSR data. BayesAss allows the estimation of symmetric or asymmetric migration rates between populations in the last two or three generations. This analysis was performed using 10,000,000 iterations, sampling every 1,000 generations and using a burn-in of 25%. The convergence of the analysis was evaluated using Tracer 1.5. The consistency of the results was evaluated by performing 10 independent runs using different seeds. RESULTS Genetic diversity and genetic-phylogeographic structure We found 26 haplotypes and 24 polymorphic sites in the control region fragment dataset. Haplotype diversities were moderate to high in all localities, being higher in the localities of the Azapa basin than in the Lluta basin; the highest haplotype diversity was found in Chapiquiña (Table 1). In the case of Socoroma, in the Lluta basin, only two haplotypes were found. The genealogy inferred by median-joining consisted of an extended haplotype network composed of two predominant haplotypes present in most localities, and multiple haplotypes of lower frequency distributed in a reticulated pattern (Fig. 2). A Socoroma private haplotype was detected, differentiated by at least six mutational steps from the rest of the network. Likewise, haplotypes similar to this Socoroma haplotype were detected in two localities of the Azapa basin (Chapiquiña and Lupica). Null alleles were detected in one of the microsatellite markers (Tchu2422); therefore, this marker was excluded, leaving a total of seven polymorphic markers which were used in the statistical analyses (Table S1). No evidence of stuttering or large allele dropout was found. No significant evidence for linkage disequilibrium was found between loci across populations. Mean expected (HE) and observed (HO) heterozygosities (microsatellites) were moderate (between 0.40 and 0.54) and of similar magnitude for most localities, except for Murmuntani, which showed lower values in both indexes. In the case of Socoroma (Lluta basin), we found a higher genetic diversity (HE) page 5 of 12 Zoological Studies 62:44 (2023) © 2023 Academia Sinica, Taiwan than in the localities of the Azapa basin (Table 1). On the other hand, a significant, negative FIS was detected in the locality of Belén, which could be explained by a sampling bias (“family bias”) in this locality. Most of the samples from Belén corresponded to larvae; therefore, only a few families would have been sampled. It has been demonstrated that the inclusion of inbred or related individuals in a sample (“family sampling”) produces a downward bias in HE (Broman 2001; Bourgain et al. 2004; De Giorgio and Rosenberg 2009; Waples 2015), which could explain the negative and significant FIS detected in Belén. Pairwise FST indices calculated both with mitochondrial (Fig. 3A) and microsatellite (Fig. 3B) markers showed that Socoroma was the most divergent locality, followed by Murmuntani and Copaquilla, showing significant values of FST in the paired comparisons (p < 0.05). The results of Geneland suggested the presence of three clusters as the most likely number of clusters (posterior density = 57.23%; average log posterior probability = -2861.1023; Fig. 4A). The most likely number of clusters was the same in all Geneland Fig. 2. Median-joining network based on the fragment of the analysed control region. Table 1. Indices of mitochondrial diversity, nuclear diversity, and inbreeding coefficients (FIS) by locality Locality n Control region SKH II n SSR HEHONAFIS Socoroma117 6 2 0.485 ± 0.079 0.00294 ± 0.00048 19 0.5316 ± 0.2027 0.5275 ± 0.2100 4.286 0.03558 Murmuntani220 6 3 0.542 ± 0.105 0.00170 ± 0.00056 22 0.3496 ± 2314 0.4026 ± 0.2643 2.286 -0.12874 Copaquilla230 8 5 0.733 ± 0.053 0.00277 ± 0.00046 32 0.4465 ± 0.2267 0.4224 ± 0.2138 3.857 0.07000 Chapiquiña326 17 9 0.883 ± 0.036 0.00494 ± 0.00064 28 0.4745 ± 0.2215 0.4434 ± 0.2223 3.714 0.08445 Belén330 8 7 0.851 ± 0.027 0.00198 ± 0.00026 30 0.4068 ± 0.2445 0.4635 ± 0.3136 3.286 -0.12240* Lupica329 16 7 0.670 ± 0.072 0.00199 ± 0.00064 32 0.4573 ± 0.2150 0.4605 ± 0.2475 4.429 0.00948 Saxamar333 13 8 0.799 ± 0.039 0.00254 ± 0.00047 36 0.4115 ± 0.2637 0.4032 ± 0.2741 3.857 0.03514 All pops. 185 24 26 0.872 ± 0.013 0.00419 ± 0.00028 199 0.4904 ± 0.1982 0.4417 ± 0.1910 6.143 0.10188* 1Socoroma River (Lluta basin); 2Seco River (Azapa basin); 3Tignamar River (Azapa basin). n = number of samples; S = number of segregating sites; K = number of haplotypes; H = haplotype diversity; Π = nucleotide diversity; HE = expected heterozygosity; HO = observed heterozygosity; NA = average number of alleles. H, Π, HE and HO are expressed as mean ± standard deviation. Significant values (p < 0.05) of FIS are indicated with an asterisk (*). page 6 of 12Zoological Studies 62:44 (2023) © 2023 Academia Sinica, Taiwan runs. The first inferred cluster was composed of the locality of Socoroma, the second cluster was composed of Murmuntani, and the third cluster included the localities of Chapiquiña, Belén, Lupica, Saxamar and Copaquilla (Fig. 4B). PCA results supported Geneland clustering; the scatter plot of the PCA (Fig. 5) showed that the Socoroma (Lluta basin) genotypes were differentiated along PC1 (7.33% of variance) from the group composed by the genotypes of the Azapa basin. Murmuntani genotypes also differed from the rest of the localities of the Azapa basin, showing high PC2 values (6.56% of variance). Gene flow Table 2 shows long-term population effective Fig. 3. Pairwise FST between localities of Telmatobius pefauri obtained using mitochondrial (A) and microsatellite (B) data. The colour scale corresponding to the values of FST is shown to the right of each matrix. Significant (Bonferroni corrected) comparisons showing p < 0.05, p < 0.01 and p < 0.001 are denoted by *, ** and ***, respectively. Fig. 4. Results of the Geneland analysis. A: Bar plot of posterior probability density according to the number of clusters; B: posterior probability maps for the delimited clusters. page 7 of 12Zoological Studies 62:44 (2023) © 2023 Academia Sinica, Taiwan sizes θ (NE*μ) and migration rates M (m*μ) obtained in MIGRATE-N (ln(Prob(D|Model)) = -1884.267708) and the number of effective migrants per generation (NE*m = θ*M). Differences were observed in the modes of the effective sizes, where θcluster1 < θcluster2 < θcluster3. Credibility intervals suggest significant differences between θcluster3 and the other two population effective sizes. Migration rate (M) and the effective number of migrants (NE*M) from cluster 2 to cluster 3 were highest. The results of the recent gene flow analysis (BayesAss) suggested a higher proportion of nonmigrants in Socoroma, Murmuntani and Saxamar Fig. 5. Scatter plot for the first two principal components obtained in the Principal Components Analysis using SSR data. Table 2. Modes and 95% credible intervals (CI) of effective size θ (NE*μ) and migration rate M (m*μ), obtained in MIGRATE-N for Telmatobius pefauri’s clusters. The number of effective migrants per generation was calculated as NE*m = θ*M Parameter Mode (95% CI) NE*m θcluster1 0.00150 (0.00000–0.00487) θcluster2 0.00223 (0.00000–0.00553) θcluster3 0.02370 (0.01660–0.03447) M 2→1 25.7 (0.0–57.3) 0.03855 M 3→1 85.0 (0.0–246.7) 0.12750 M 1→2 0.3 (0.0–435.3) 0.00067 M 3→2 0.3 (0.0–422.0) 0.00067 M 1→3 0.3 (0.0–302.0) 0.00711 M 2→3 457.7 (35.3–492.7) 10.84749 page 8 of 12Zoological Studies 62:44 (2023) © 2023 Academia Sinica, Taiwan (Table 3). Higher mean migration rates (with standard deviations below 0.05) were detected from Saxamar to Chapiquiña, Belén and Lupica, all belonging to the Tignamar River. Migration from Belén to Copaquilla (Seco River) was also detected, but showed a high standard deviation (> 0.1). DISCUSSION Regarding the initial predictions, the results of the analyses of genetic structure and gene flow suggest that populations of T. pefauri belonging to the same drainage (Tignamar River) were less differentiated (nonsignificative FST values in all except one comparison) and exhibited higher gene flow levels than their conspecifics inhabiting the other drainage in the same basin (Seco River), and those belonging to the other basin (Lluta basin/Socoroma). In addition, populations of T. pefauri located in different basins were more differentiated (higher FST values) than their conspecifics from different drainages, but the same basin. This hierarchical pattern of genetic structure more closely matches the predictions of SHM than the “random” pattern of differentiation predicted by the DVM, which is characterized by a complete lack of spatial genetic structure explained by landscape structure, such as catchment boundaries or even overall geographic distance among populations, and the lack of gene flow between populations (Hughes et al. 2009). Nevertheless, high differentiation was found between Murmuntani and Copaquilla, which belong to the same drainage, the Seco River. Indeed, Murmuntani constituted a different cluster within the Azapa basin. Considering that the highest levels of differentiation were found between populations from different basins, the Murmuntani cluster could be explained by the geographical isolation and peripheric position of this locality in relation to the other T. pefauri’s populations within the Azapa basin. Murmuntani belongs to the Seco River drainage, which receives little precipitation compared to the Tignamar River drainage (Niemeyer 1982). Therefore, higher levels of fragmentation in the Seco River could explain the differentiation of the Murmuntani cluster. The other known population of T. pefauri present in the Seco River drainage, Copaquilla, is located near the confluence of the Seco and Tignamar Rivers and showed more similarity to the populations inhabiting the Tignamar River drainage than to Murmuntani, supporting the idea of geographical isolation in Murmuntani. The contemporary gene flow detected between the tributaries of the Tignamar River would be asymmetric, where migration would have occurred from Saxamar to the rest of the localities, congruent with the geographic arrangement of the tributaries (i.e., Saxamar is the southernmost locality of the distribution) and the northern direction of the water flow in the Tignamar River. However, sampled populations of the Tignamar River drainage are located in separate tributaries, not the main river. If individuals migrate from Saxamar to Table 3. Migration rates (recent) between localities inferred with BayesAss. Posterior mean migration rates and the standard deviation of the marginal posterior distribution for each estimate (in parentheses) are shown. The localities from which each individual was sampled are listed in the rows, while the localities from which they migrated are listed in the columns. Values along the diagonal are the proportions of individuals derived from the source populations in each generation. Migration rates > 0.10 are shown in bold Drainage Socoroma river Seco river Seco river Tignamar river Tignamar river Tignamar river Tignamar river Locality 1234567 1 – Socoroma 0.8652 (0.0365) 0.0137 (0.0132) 0.0239 (0.0211) 0.0154 (0.0151) 0.0224 (0.0200) 0.0443 (0.0287) 0.0150 (0.0143) 2 – Murmuntani 0.0115 (0.0112) 0.9031 (0.0342) 0.0126 (0.0123) 0.0119 (0.0116) 0.0170 (0.0165) 0.0128 (0.0126) 0.0310 (0.0245) 3 – Copaquilla 0.0089 (0.0087) 0.0105 (0.0103) 0.7348 (0.0857) 0.0137 (0.0133) 0.1795 (0.1042) 0.0133 (0.0131) 0.0392 (0.0369) 4 – Chapiquiña 0.0094 (0.0090) 0.0341 (0.0258) 0.0151 (0.0164) 0.6793 (0.0126) 0.0718 (0.0470) 0.0151 (0.0149) 0.1752 (0.0501) 5 – Belén 0.0091 (0.0088) 0.0269 (0.0227) 0.0157 (0.0160) 0.0115 (0.0114) 0.6868 (0.0194) 0.0163 (0.0160) 0.2337 (0.0363) 6 – Lupica 0.0089 (0.0087) 0.0199 (0.0176) 0.0179 (0.0179) 0.0149 (0.0146) 0.0489 (0.0377) 0.6830 (0.0156) 0.2064 (0.0430) 7 – Saxamar 0.0079 (0.0076) 0.0212 (0.0185) 0.0212 (0.0199) 0.0124 (0.0117) 0.0463 (0.0363) 0.0239 (0.0222) 0.8671 (0.0455) page 9 of 12Zoological Studies 62:44 (2023)