Old divergence and restricted gene flow between torrent duck (Merganetta armata) subspecies in the Central and Southern Andes

Abstract Aim To investigate the structure and rate of gene flow among populations of habitat‐specialized species to understand the ecological and evolutionary processes underpinning their population dynamics and historical demography, including speciation and extinction. Location Peruvian and Argentine Andes. Taxon Two subspecies of torrent duck (Merganetta armata). Methods We sampled 156 individuals in Peru (M. a. leucogenis; Chillón River, n = 57 and Pachachaca River, n = 49) and Argentina (M. a. armata; Arroyo Grande River, n = 33 and Malargüe River, n = 17), and sequenced the mitochondrial DNA (mtDNA) control region to conduct coarse and fine‐scale demographic analyses of population structure. Additionally, to test for differences between subspecies, and across genetic markers with distinct inheritance patterns, a subset of individuals (Peru, n = 10 and Argentina, n = 9) was subjected to partial genome resequencing, obtaining 4,027 autosomal and 189 Z‐linked double‐digest restriction‐associated DNA sequences. Results Haplotype and nucleotide diversities were higher in Peru than Argentina across all markers. Peruvian and Argentine subspecies showed concordant species‐level differences (ΦST mtDNA = 0.82; ΦST autosomal = 0.30; ΦST Z chromosome = 0.45), including no shared mtDNA haplotypes. Demographic parameters estimated for mtDNA using IM and IMa2 analyses, and for autosomal markers using ∂a∂i (isolation‐with‐migration model), supported an old divergence (mtDNA = 600,000 years before present (ybp), 95% HPD range = 1.2 Mya to 200,000 ybp; and autosomal ∂a∂i = 782,490 ybp), between the two subspecies, characteristic of deeply diverged lineages. The populations were well‐differentiated in Argentina but moderately differentiated in Peru, with low unidirectional gene flow in each country. Main conclusions We suggest that the South American Arid Diagonal was preexisting and remains a current phylogeographic barrier between the ranges of the two torrent duck subspecies, and the adult territoriality and breeding site fidelity to the rivers define their population structure.

In the rivers occupied by torrent ducks, the adults show a high breeding site fidelity (Alza et al., 2017;Cardona & Kattan, 2010) to territories that reportedly range from 0.7 to 2 km during the breeding season (Cerón & Trejo, 2009;Colina, 2010;Eldridge, 1986;Johnsgard, 2010;Moffet, 1970;Naranjo & Ávila, 2003;Ubeda, Cerón, & Trejo, 2007). However, immature and adult individuals have been reported moving away from rivers likely searching for new habitats (i.e., four opportunistic records of male-biased dispersal ;Cerón & Capllonch, 2016). Thus given that torrent duck pairs are territorial with a high rate of river-specific site fidelity, we expect torrent ducks to be more genetically similar among the nearest rivers but isolated by distance among rivers depending on the dispersal capacity of the ducks (Hanski & Gilpin, 1991;Harrison & Hastings, 1996;Nielsen & Slatkin, 2013;Wright, 1943). Alternatively, different riverine subpopulations may follow a metapopulation model (Hanski & Gaggiotti, 2004;Levins, 1968), in which populations or subpopulations isolated including no shared mtDNA haplotypes. Demographic parameters estimated for mtDNA using IM and IMa2 analyses, and for autosomal markers using ∂a∂i (isolation-with-migration model), supported an old divergence (mtDNA = 600,000 years before present (ybp), 95% HPD range = 1.2 Mya to 200,000 ybp; and autosomal ∂a∂i = 782,490 ybp), between the two subspecies, characteristic of deeply diverged lineages. The populations were well-differentiated in Argentina but moderately differentiated in Peru, with low unidirectional gene flow in each country.

Main conclusions:
We suggest that the South American Arid Diagonal was preexisting and remains a current phylogeographic barrier between the ranges of the two torrent duck subspecies, and the adult territoriality and breeding site fidelity to the rivers define their population structure.

K E Y W O R D S
Andes, gene flow, genetic diversity, Merganetta armata, population structure, time since divergence inside watersheds are still connected by limited gene flow, that would play a key role in recolonization (sporadic dispersal) after the extinction of a subpopulation.
Here, we reconstruct population structure for mitochondrial DNA (mtDNA) and thousands of nuclear markers isolated using partial genome resequencing (i.e., double-digest restriction site-associated DNA sequencing (ddRAD-Seq)), attained from altitudinal sampling transects in Peru and Argentina. Specifically, we estimated the time since divergence and gene flow between the two subspecies, and also characterized population structure and gene flow between riverine populations within subspecies. Furthermore, we calculated genetic diversity across the different populations to infer whether specific riverine populations have experienced demographic processes such as population expansion or contraction. Because the spatial linear configuration of the rivers restricts suitable habitat and constrains the population, we also examined evidence for altitudinal clines in mtDNA haplotype frequency.
Finally, this genetic characterization across altitudinal gradients permits us to explore whether genetic variation is associated with recently described morphological and physiological differences of torrent ducks at high elevations, including changes in body size (Gutiérrez-Pinto et al., 2014), function of key enzymes (Dawson et al., 2016), changes in insulative properties of feathers (Cheek, Alza, & McCracken, 2018), and changes in the hypoxic ventilatory response (Ivy et al., 2019).

| Sampling and DNA extraction
Blood samples were collected from 156 captured and released torrent ducks, during the 2010-2011 dry seasons, from two rivers in the western and eastern slopes of the Central Andes of Peru (Chillón River, n = 57 and Pachachaca River, n = 49; Table 1 and Figure 1b,c), and two rivers in the eastern slope of the southern Andes of Argentina (Arroyo Grande River, n = 32 and Malargüe River, n = 18; Table 1; Figure 1b,c). Each river was surveyed along an altitudinal transect (maximum range, 900-4,000 m) using an active mist net method by which torrent ducks were driven to nets (Alza et al., 2017). Blood samples (3 ml of whole blood per individual) were stored in liquid N 2 in the field before being placed in long-term storage at −80°C. Blood samples are archived at the University of Alaska Museum (Fairbanks, Alaska), Museo Argentino de Ciencias Naturales "Bernadino Rivadavia" (Buenos Aires, Argentina), and CORBIDI (Lima, Peru). DNA was extracted using a DNeasy Blood & Tissue kit and following the manufacturer's protocols (Qiagen). Extractions were quantified using a NanoDrop 2000 Spectrophotometer (Thermo Fisher Scientific Inc.) to ensure a minimum concentration of 20 ng/µl.

| ddRAD-Seq library preparation
While mtDNA was obtained across all samples, for more extensive genomic analysis of the two subspecies, a subset from Chillón River in Peru (n = 10) and Malargüe River in Argentina (n = 9) were subjected to partial genome resequencing via ddRAD-Seq. Sample preparation for ddRAD-Seq followed the double-digest protocol outlined in DaCosta and Sorenson (2014). In short, ~1 μg of genomic DNA was digested using 10 U of each SbfI and EcoRI restriction enzymes. Adapters containing sequences compatible for Illumina TruSeq reagents and barcodes for de-multiplexing reads were ligated to the sticky ends generated by the restriction enzymes.
The adapter-ligated DNA fragments were then size-selected for 300-450 bp using gel electrophoresis (2% low-melt agarose) and purified using a MinElute Gel Extraction Kit (Qiagen). Size-selected fragments were then PCR amplified with Phusion High-Fidelity DNA Polymerase (Thermo Scientific), and the amplified products were cleaned using AMPure XP magnetic beads (Beckman Coulter, Inc.

| Bioinformatics of ddRAD-Seq data
Raw Illumina reads were processed using a computational pipeline described by DaCosta and Sorenson (2014; http://github.com/BU-RAD-seq/ddRAD-seq-Pipeline). First, reads were assigned to individual samples based on barcode sequences using the ddRADparser.
py script. Reads per sample were then collapsed into identical clusters using the CondenseSequences.py script with low-quality reads (i.e., sequences that failed to cluster with any other reads (-id setting of 0.90) and with an average per-base Phred score <20) filtered out with the FilterSequences.py script. Condensed and filtered reads from all samples were then concatenated and clustered with anid setting of 0.85 in ucluST. MuScle v.3 (Edgar, 2004) was used to align and cluster reads, and samples within each aligned cluster were genotyped using the RADGenotypes.py script. Homozygotes and heterozygotes were identified based on thresholds outlined in DaCosta and Sorenson (2014;Lavretsky et al., 2015), with individual genotypes falling into three categories: "missing" (no data), "good" (unambiguously genotyped), and "flagged" (recovered heterozygous genotype, but with haplotype counts outside of acceptable thresholds or with >2 alleles detected). Loci with <20% missing genotypes and ≤6 flagged genotypes were retained for downstream analyses.
Moreover, alignments with end gaps due to indels, ≥2 polymorphisms in the last five base pairs, and/or a polymorphism in the increased the total number of retained markers by ~7%, while reducing bias resulting from discarding loci with indels or high levels of polymorphisms. Finally, output files for downstream analyses were created with custom python scripts ) that incorporated sequencing depth. To limit any biases due to sequencing error and/or low depth, alleles were called "missing" unless they met our thresholds of a minimum of 5X coverage for homozygotes, and thus at least 10X coverage for heterozygotes to call both alleles with >99% predicted sequence accuracy.
Finally, autosomal and Z chromosome-linked loci were identified as described in Lavretsky et al. (2015), with assignments based on differences in sequencing depth and homozygosity between males and females. Because females have only one Z chromosome, Zlinked markers in females are expected to appear homozygous and be recovered at about half the sequencing depth of males.

| Summary statistics
For mtDNA, we calculated summary statistics between rivers including haplotype diversity, nucleotide diversity (π/site), Tajima Analysis of molecular variance (AMOVA) was run to examine genetic differentiation within and among the four rivers sampled from the two countries. Population pairwise Φ ST was estimated under the Tamura and Nei (1993) model of nucleotide substitution. For ddRAD-Seq autosomal and Z chromosome-linked markers, pairwise Φ ST (i.e., "nuc.F_ST") and nucleotide diversity (i.e., "nuc.div.within") estimates were calculated in the R program "PopGenome" (Pfeifer, Wittelsbürger, Ramos-Onsins, & Lercher, 2014); indel positions were excluded from analyses.

| Population structure
For mtDNA, the population structure was visualized by reconstructing a haplotype network using TcS v.1.21 (Clement, Posada, & Crandall, 2000). TcS illustrates all connections that have a 95% probability of being the most parsimonious. Networks are more appropriate for intraspecific gene genealogies than rooted tree algorithms because population genealogies are often multifurcated, descendant genes coexist with persistent ancestral sequences and in the case of nuclear DNA recombination events produce reticulate relationships, where traditional phylogenetic trees treat all sequences as terminal taxa (Posada & Crandall, 2001).
Analysis of population subdivision for ddRAD-Seq data was done using two methods. First, we used a principal component analysis (PCA) as implemented in the adegenet R program (i.e., "dudi.pca"; Dray & Dufour, 2007; also see Jombart, 2008).

| Historical population demography of nuclear DNA
We further estimated rates and directionality of gene flow from the ddRAD-Seq data with the program ∂a∂i (Gutenkunst, Hernandez, Williamson, & Bustamante, 2009. ∂a∂i implements an efficient diffusion approximation-based approach to test empirical data against specified evolutionary models (e.g., isolation-with-migration). Using ∂a∂i, a site frequency spectrum was derived from all biallelic RAD-Seq autosomal SNPs. Loci were concatenated, and SNPs extracted and formatted for ∂a∂i using the custom python script nex_mo.py . Because we lacked an outgroup, site frequency spectrum data were folded, with only minor alleles considered in the frequency spectrum. Variants observed in zero or in all samples were ignored (masked), as described by Gutenkunst, Hernandez, Williamson, and Bustamante (2010). Finally, for ∂a∂i to accommodate missing data and differences in sample sizes be- To convert the parameter estimates from ∂a∂i to biologically informative values, we estimated generation time (G) and mutation rates (μ, per locus). First, generation time (G) was calculated as G = α + (s/(1 − s)), where α is the age of maturity and s is the expected adult survival rate (Saether et al., 2005). Although sexually active by the first generation, most duck species reach sexual maturity in their second year (α = 2) with an average adult survival rate of 0.37 estimated for torrent ducks from mark-recaptured data of the Chillón River (pers. obs.). Together, we estimated a generation time to be 2.5 years. Next, to obtain a mutation rate for nuclear genes, we multiplied a rate of 1.2 × 10 -9 substitutions/site/year, previously calculated for nuclear genes in other ducks (Peters, Zhuravlev, Fefelov, Humphries, & Omland, 2008) by generation time to attain a rate of 3 × 10 -9 substitutions site −1 generation −1 (s s −1 g −1 ). A final mutation rate was calculated as the product of the above mutation rate and the total number of base pairs.

| Genetic diversity & population structure-mtDNA
Nucleotide (π/site) and haplotype diversities of the mtDNA CR differed among riverine populations and countries (Table 2).
M. a. leucogenis in Peru had higher nucleotide diversity that was five times greater than M. a. armata in Argentina, although the TA B L E 2 Nucleotide and haplotype diversity for mtDNA control region (Chillón River, n = 57; Pachachaca River, n = 49; Arroyo Grande River, n = 33 and Malargüe River, n = 17), autosomal and Z-linked markers (Peru-Chillón, n = 10 and Argentina-Malargüe, n = 9) from four populations and two subspecies of torrent duck (Merganetta armata) in the Andes of Peru and Argentina  Table 2, Table S1). In Argentina, we once again found a single mtDNA haplotype shared between torrent ducks from the Arroyo Grande and Malargüe rivers, with the Arroyo Grande River having a higher number of private haplotypes as compared to those from the Malargüe River ( Figure 1a; Table 2,   Table S1). Despite recovered structure among torrent duck populations (Table 3) Figure 2a). A single biallelic SNP was randomly chosen across loci, with a total of 2,356 and 49 biallelic autosomal and Z-linked SNPs, respectively, used for ADMIXTURE analyses.

| Genetic diversity & population structure-ddRAD-seq loci
The optimal K was two for both autosomal and Z-linked markers, and ADMIXTURE results largely distinguished between the two subspecies with 99% of assignment probabilities and were concordant with the PCA results, plotted with the first two principal components, in which the two subspecies were differentiated in two discrete clusters, whereas individuals from Peru are much more dispersed along PC1 than the individuals from Argentina (Figure 2b; Figure S1). No additional structural resolution was attained when analyzing higher values of K.

| Historical population demography
The effective population size parameters, Θ, scaled to the neutral mutation rate in IM and IMa2 for mtDNA were higher for the rivers in Peru than Argentina (Figures 3b and 4b). The largest Θ was for the Pachachaca River with a peak at 19.3 (95% HPD = 9.9-35.0;  (Table 3) were not significantly different from zero (higher p-values). Therefore, we did not observe a significant excess of rare nucleotide site variants compared to what would be expected under a neutral model of evolution, and each riverine population may be evolving as per mutation-drift equilibrium.
In contrast, Tajima's D values were negative for autosomal and

Z-linked markers for both subspecies populations in Peru and
Argentina, although both markers had not statistically significant p-values (Table 3).
Estimating gene flow, the four-population IMa2 analysis suggests that the number of migrants (m) between riverine populations within each country is effectively low and unidirectional ( Figure 3c).
Importantly, we were unable to reject the hypothesis of no gene flow. For time since divergence, the posterior probability distribution of t (scaled time since divergence) with the IMa2 analysis between the two torrent duck subspecies peaked at 17.32 (95% HPD = 9.3-24.5; Figure 2d). This value converted to time in years (Peters et al., 2005)

| D ISCUSS I ON
We provide the most comprehensive molecular analysis of torrent ducks to date and clearly identify limited or no gene flow and strong differentiation across the mitochondrial and nuclear genomes of torrent ducks found in the rivers of Peru versus Argentina. These results support the current taxonomic subspecies designation of Peruvian and Argentine torrent ducks. Moreover, we are able to conclude that these two subspecies probably originated during the early-to mid-Pleistocene as a result of an old divergent nonvicariant event subsequent to the origin of the South American Arid Diagonal phylogeographic barrier (Burniard, 1982). Finally, we also recovered strong population structure within each country that is probably explained by both a strong territoriality and breeding site fidelity in adults, and the geographic isolation by distance between riverine subpopulations.

| Genetic diversity and small effective population size
In general, diversity estimates across mtDNA haplotypes for the two torrent ducks subspecies were higher or similar to those found in other Andean duck species (  Hourlay et al., 2008)). In fact, torrent ducks from Peru had values of mtDNA diversity more similar to those found in other ducks with Ne several orders of magnitude larger Ne (e.g., mallards; Lavretsky et al., 2015). Therefore, because mutation rates are likely to be similar across duck species, we can hypothesize that torrent duck mtDNA genetic diversity probably is influenced by the structure of their habitat leading to extensive branched and subdivided, structured riverine habitat, such that each river may be isolated by distance and genetically distinct from other watersheds. In contrast to mtDNA, diversity estimates from nuclear DNA for the two torrent duck subspecies were lower compared to other waterfowl species (Lavretsky et al., 2015;Lavretsky et al., 2016;Peters et al., 2016;Wilson et al., 2012). A probable explanation for the disparity in diversity estimates of mitochondrial and nuclear genes within torrent ducks could be a combined effect resulting from high mtDNA mutation rate (four orders of magnitude), population subdivision due to the branched structure riverine habitat, the adult territorial behavior, and natural selection associated to high elevation (e.g., low temperatures and hypoxia) affecting nuclear DNA.
Finally, between countries (Argentina vs. Peru) and within countries (e.g., Chillón River vs. Pachachaca River), we found that populations inhabiting shorter rivers compared to populations inhabiting longer or branched rivers systems had the lowest levels of molecular diversity (Figure 1; Tables 1 and 2). This low genetic diversity associated with river length is more evident in the southern latitudinal rivers sampled in Argentina, which may have experienced founder events following glaciation periods. In general, these rivers are shorter than the Peruvian rivers and more likely maintaining small Ne that might be more susceptible to loss of diversity due to genetic drift (Table 1; Figure 3b). Moreover, current and historical changes in the river water regimes due to extreme drought and dramatic climate changes (Rivera, Penalba, Villalba, & Araneo, 2017) can produce bottlenecks, or riverine population extinctions. These  Table 3).

| Deep genetic divergence between subspecies
We found mtDNA monophyly and nonsignificant estimates of gene flow between the two torrent duck subspecies, M. a. leucogenis and M. a. armata (Figures 1a and 4a,c; Table 3; also for nuclear DNA Figure 2b, Figure S1). The reduced gene flow between the two torrent ducks subspecies is most likely due to their habitat specialization and lack of suitable intermediate habitat, particularly where the Andes show the highest mountain ranges and widest extension that harbors a complex riverine spatial structure (Capitonio et al., 2011;Gonzalez & Pfiffner, 2011). In contrast, two subspecies of speckled teal (Anas flavisrostris) with similar distribution range and divergence time show relatively less genomic divergence and asymmetric gene flow (Graham et al., 2017). Long-standing and recurrent gene flow in speckled teal might be facilitated by the lack of habitat specialization (i.e., lakes, pond, and also rivers), and a larger number of habitats present in the lowlands due to a broader flat topography than in the highlands in the Southern Andes. Therefore, when and where geographic barriers isolate these torrent duck subspecies is likely tightly connected with the orogeny of the Andes, and the formation and persistence of torrential riverine systems (Fjeldså, 1985).
Finally, we estimate that the two torrent ducks likely diverged in the early-to mid-Pleistocene (mtDNA divergence time = 600,000 ybp, 95% HPD range = 1.2 Mya to 200,000 ybp, Figure 3d; and nuclear divergence time = 782,490 ybp). Compared to other species found in the Andes of Argentina and Peru, divergence estimates for torrent ducks are some of the oldest, although comparable to those estimated between speckled teal (Bulgarella et al., 2012;Graham et al., 2017;McCracken et al., 2009;Wilson et al., 2012). Therefore, torrent ducks can be considered an older Andean resident probably tightly associated to the formation of river systems. However, divergence of these two subspecies likely followed the establishment of all major river systems as the Andes started to rise during the Late Cretaceous (~60-80 Mya) and were at their final heights by the Early Pliocene (~3-4 Mya) (Capitonio et al., 2011;Clapperton, 1983;Garzione et al., 2008;Hartley, 2003), which are well before estimated divergence times.
Moreover, given that these historical climate conditions of the Arid Diagonal likely reduced the river habitat suitability long before the estimated time since divergence between these two torrent duck subspecies, we hypothesize that colonization through dispersal and subsequent isolation likely best explains their divergence.
Although rare, torrent ducks are known to make long distance movements (Cerón & Capllonch, 2016) that would permit them to traverse across large landscapes like the Dry Diagonal. In general, we hypothesize that it is very likely that these two subspecies were probably originated by an old divergence through multiple founder events, with individuals dispersing from the Central Andes to the Southern Andes, because the large Ne and genetic diversity recovered in Peru compared to Argentina, during the early-to mid-Pleistocene.

| Population structure among rivers, "extended family" within rivers and metapopulation
Based on mitochondrial DNA, we estimated low or zero gene flow between rivers within each country, suggesting well-differentiated and moderately structured river populations in Peru and Argentina, respectively. These results indicate that each river system, inside of a watershed, contains a particular group of genetically more similar individuals. We hypothesize that these results are likely due to three nonmutually exclusive hypotheses: (1) Strong territoriality and breeding site fidelity in adults increase the similarity of the individuals within a river, (2) geographic (e.g., watershed boundary, high mountain ranges) or climatic barriers (e.g., large deserts) isolate rivers preventing the dispersal of the individuals, and driving an independent evolution of each riverine population by genetic drift, and (3) individual dispersal capacity (i.e., frequency and distance) describes an isolation-by-distance pattern that in short distance dispersal resembles a metapopulation model. For example, the year-round territorial behavior, the long-term pair bonds, and the strong site fidelity in adult torrent ducks, as mentioned in the first hypothesis, can prevent gene flow among rivers, but also contribute to the similarity of the individuals of a growing population in vacant or improved habitat, and thus increase the relatedness within each riverine population and develop an "extended family" (Alza et al., 2017;Cardona & Kattan, 2010;Eldridge, 1986;Hartl & Clark, 2007;Pernollet, Estades, & Pavez, 2012). Additionally, as referred in the second hypothesis, torrential river systems inside watersheds are analogous to islands (islands on mountain slopes) that isolate and sustain populations and communities (Black, 1997;Naiman, Magnuson, McKnight, Stanford, & Karr, 1995;Omernik & Bailey, 1997;Sullivan, Watzin, & Keeton, 2007), similarly observed in lakes (Barbour & Brown, 1974), marshes (Brown & Dinsmore, 1988), caves (Culver, Holsinger, & Baroody, 1973), mountaintops (Nores, 1995), or woodlots (Holland, 1978). Thus, as new rivers are colonized (e.g., by a small group of individuals or for a long term), these isolated populations may behave independently of the other or larger population, growing as distinct populations and leading to increased population structure among rivers (Mayr, 1963;Phillimore & Owens, 2006). Finally, it is possible that different watershed (sub)-populations may be genetically more connected or completely differentiated depending on the distance among them and the dispersal capacity of the ducks. Among torrent duck riverine (sub)-populations it appears that dispersal is most likely to occur at close range (Cerón & Capllonch, 2016;Paradis, Baillie, Sutherland, & Gregory, 1998), which could also reduce their differentiation and resemble a metapopulation model (Hanski & Gaggiotti, 2004;Levins, 1968). Similar population structure patterns have been described in the blue duck (Hymenolaimus malacorhynchos) in New Zealand with a regional genetic differentiation (region around 100 km) associated with isolation-by-distance pattern (Grosser et al., 2016). To test these hypotheses, future work would benefit from measures of relatedness on increased sampling of individuals across rivers, including multiple adjacent river systems stratified by elevation. We note that complete haplotype admixture (mtDNA) along each study river system suggests that torrent ducks are not segregated by elevation (considering a 2,500 m divide, Figure 1b).
Previous studies have reported seasonal altitudinal movements of torrent ducks, and these movements can be influenced by reproductive season, food availability, variation in the streamflow, and winter conditions, especially in the Southern Andes (Johnsgard, 1966;Johnson, 1963;Pernollet et al., 2012;Ramírez, Botero, & Kattan, 2014). Also during our banding activities in the Chillón River, we recovered two individuals that had moved at least 4 km along the river.

| Conservation relevance
Torrent ducks are vulnerable to systematic pressure and stochastic perturbations because of their specialization to whitewater streams, low population densities, patchy distribution, low reproductive potential, and reliance on water quality and aquatic insect populations (Alvarez, Astie, Debandi, & Scheibler, 2014;Carboneras, 1992;Shaffer, 1981;Zwick, 1992). Indeed, populations living in rivers and streams habitats are susceptible to extinction as a result of the high vulnerability that freshwater systems have pollution, siltation, and anthropogenic disturbance. Thus, torrent ducks had suffered local extinction and decreasing population trends of many river populations (Callaghan, 1997;Cerón & Trejo, 2012;Pernollet, Pavez, & Estades, 2013;Weller, 1975). In particular, the recognized causes of their riverine extinction are as follows: predation on offspring by exotic invasive species (mink; Cerón & Trejo, 2012), competition for aquatic insect food sources (trout; Eldridge, 1986), loss of territories and habitats due to urbanization and contamination of rivers in Peru and Colombia (J. M. Barbarán, personal communication; Cardona & Kattan, 2010), management of the watershed to provide drinking water, irrigation and hydropower production for mining activities and human communities (Pernollet et al., 2013), and hunting without proper regulation (pers. obs.). Whereby, the three subspecies regionally had been classified under threatened categories (Vulnerable and Endangered) (Cerón & Trejo, 2012;Green, 1996), even though the overall species is still classified as a "Least Concern" (BirdLife International, 2016). Our study also emphasizes the necessity to protect and manage the subspecies and riverine populations as separate conservation units due to the strong genetic structure and very low gene flow. Therefore, the establishment of monitoring programs of the subspecies and multiple populations is a priority to evaluate their metapopulation structure and understand the actual status of the species (Phillimore & Owens, 2006), even more, under the threat of deglaciation and regional drought conditions associated to climate change in the Andes. In general, the torrent duck provides an excellent opportunity to examine the patterns and the mechanisms related to strong population structure among small and isolated populations, and these results can be contrasted with other riverine ducks, such as the African black duck (Anas sparsa), and Salvadori's duck (Salvadorina [Anas] waigiuensis) in New Guinea, or habitat specialist species.