Neutral and functionally important genes shed light on phylogeography and the history of high‐altitude colonization in a widespread New World duck

Abstract Phylogeographic studies often infer historical demographic processes underlying species distributions based on patterns of neutral genetic variation, but spatial variation in functionally important genes can provide additional insights about biogeographic history allowing for inferences about the potential role of adaptation in geographic range evolution. Integrating data from neutral markers and genes involved in oxygen (O2)‐transport physiology, we test historical hypotheses about colonization and gene flow across low‐ and high‐altitude regions in the Ruddy Duck (Oxyura jamaicensis), a widely distributed species in the New World. Using multilocus analyses that for the first time include populations from the Colombian Andes, we also examined the hypothesis that Ruddy Duck populations from northern South America are of hybrid origin. We found that neutral and functional genes appear to have moved into the Colombian Andes from both North America and southern South America, and that high‐altitude Colombian populations do not exhibit evidence of adaptation to hypoxia in hemoglobin genes. Therefore, the biogeographic history of Ruddy Ducks is likely more complex than previously inferred. Our new data raise questions about the hypothesis that adaptation via natural selection to high‐altitude conditions through amino acid replacements in the hemoglobin protein allowed Ruddy Ducks to disperse south along the high Andes into southern South America. The existence of shared genetic variation with populations from both North America and southern South America as well as private alleles suggests that the Colombian population of Ruddy Ducks may be of old hybrid origin. This study illustrates the breadth of inferences one can make by combining data from nuclear and functionally important loci in phylogeography, and underscores the importance of complete range‐wide sampling to study species history in complex landscapes.


| INTRODUC TI ON
Phylogeographic analyses aimed at understanding historical processes underlying species distributions typically focus on associations between geography and the distribution of neutral genetic variation (Avise, 2000;Hickerson et al., 2010;Knowles, 2004).
However, genes of functional importance can also be employed to investigate the history of species, particularly when populations are exposed to distinct environments with contrasting selective pressures (Deagle, Jones, Absher, Kingsley, & Reimchen, 2013;Feldman, Brodie, & Pfrender, 2009;Hoekstra, Drumm, & Nachman, 2004;Savage & Zamudio, 2016;Sork et al., 2016). Analyses of functionally important genes allow one to address questions about evolutionary processes such as adaptation and natural selection, complementing inferences of demographic processes based on patterns of geographic variation in neutral genes (Zamudio, Bell, & Mason, 2016). Here, we integrate analyses of neutral markers and proteincoding genes potentially under selection to test hypotheses about the biogeographic and evolutionary history of a widespread species of duck with a broad latitudinal and elevational distribution in the New World.
A long-standing question in biogeography is that of the evolutionary origin of high-elevation organisms. In the New World, in particular, researchers have long been interested in understanding the origins of biotas occurring in the high Andes of South America, with a leading hypothesis being that tropical high-elevation species are derived from low-elevation ancestors occupying areas with similarly cool climates in the temperate zone (Chapman, 1917;Vuilleumier, 1986b). Examples of plant and animal lineages colonizing the high Andes from temperate regions and subsequently diversifying in tropical montane areas are consistent with the hypothesis that such lineages were preadapted to occur in cool high-elevation environments and that this allowed them to occupy newly formed habitats resulting from the Andean uplift (Hughes & Eastwood, 2006;Sanín et al., 2009). However, because life at high elevations requires not only dealing with low temperatures but also coping with hypoxia, colonization of montane environments by species from lowelevation areas may have also entailed adaptive evolution in traits related to oxygen (O 2 )-transport physiology.
High-elevation environments are of particular interest for research at the interface of biogeography, adaptation, and physiology because hypoxia increases markedly along elevational gradients, likely imposing strong selective pressures and setting limits to species ranges. In high-elevation regions above 4,000 m, such as the high Andes, the partial pressure of oxygen is approximately 60% of that at sea level (Beall, 2006;Hopkins & Powell, 2001). Under such conditions, oxygen transport to tissues is reduced, altering metabolism and inhibiting physical activity (Hopkins & Powell, 2001;Storz & Moriyama, 2008). A variety of physiological compensatory mechanisms have evolved in high-altitude organisms across the oxygen transport cascade, examples of which include both genetic adaptations and environmentally induced plasticity (Storz, 2010). One of the best-studied examples of genetically based adaptation involves the blood O 2 transport protein hemoglobin. Many organisms native to high altitude exhibit hemoglobin with a high affinity for oxygen (Weber, 2007); in birds, increased hemoglobin-oxygen affinity has evolved convergently across multiple lineages occurring at high elevations (Natarajan et al., 2016). Studies of hemoglobin have also revealed environmentally induced phenotypic plasticity in response to pressures associated with high elevations. Unlike many high-altitude endemics, which often maintain relatively low hemoglobin concentrations, species native to low altitude typically increase their hemoglobin concentration and associated hematocrit (i.e., the volume of packed red blood cells) upon ascending to high altitude, thereby increasing their total blood oxygen carry capacity (Lui et al., 2015;Monge & Leon-Velarde, 1991;Tufts et al., 2013). Although there is evidence indicating that such response may be maladaptive over the long term due to increased blood viscosity (Storz & Moriyama, 2008;Storz, Scott, & Cheviron, 2010), over the short term, this plasticity is an essential component of the acclimatization response of many low-altitude organisms.
Analyses of geographic variation in genes involved in O 2 transport have provided important insights about the history of occupation of high-elevation environments by birds. For example, estimates of gene flow between lowland and highland populations of Crested Ducks (Lophonetta specularioides) obtained based on hemoglobin alleles are much lower than those inferred from neutral markers; in light of biochemical properties of amino acid substitutions differing across elevations, such differences in hemoglobin gene flow are consistent with a role for selection structuring genetic and phenotypic variation (Bulgarella et al., 2012). Likewise, hemoglobin variants with increased affinity for oxygen occur at a greater frequency in highelevation populations of House Wrens (Troglodytes aedon) relative to those from lowland areas, and F st across elevations for a functionally important site in the βA-globin gene is greater than F st values for most coding SNPs across the genome (Galen et al., 2015). Because the high-elevation βA-globin genotype is derived, colonization of high-elevation environments by House Wrens presumably involved adaptive evolution in protein function.
A case in which adaptive evolution of hemoglobins has been hypothesized to play a role in the evolutionary history of a species and facilitated its spread across landscapes is that of the Ruddy Duck (Anatidae, Oxyura jamaicensis), a diving duck found throughout the New World from southern Canada to Tierra del Fuego (Figure 1).
Three subspecies are recognized based on facial plumage of males (Adams & Slavid, 1984;Fjeldså, 1986;Siegfried, 1976 Furthermore, based on geographic variation in substitutions in the βA globin gene, it was proposed that, upon expanding its range, the species first acclimatized or adapted to highlands of the Andes (Muñoz-Fuentes et al., 2013). A substitution putatively increasing the affinity of hemoglobin for oxygen (McCracken, Barger, Bulgarella, Johnson, Sonsthagen, et al., 2009) may have allowed Ruddy Ducks to establish in the Northern Andes and to spread southward along highland environments; as the species colonized precordilleran steppe habitats in southern South America, its hemoglobin seemingly adapted again to lowland conditions (Muñoz-Fuentes et al., 2013). These inferences, however, were made based on molecular data sets that did not include population-level sampling in the Northern Andes and were not validated by functional analyses. In fact, a subsequent analysis revealed that substitutions occurring in high-elevation Ruddy Ducks have no discernible effects on the oxygen affinity of hemoglobin . Although this analysis did not rule out the possibility that genetic differences between lowland and highland Ruddy Ducks may affect other structural or functional properties of hemoglobins , it implies that the previously Here, we revisit the phylogeography of Ruddy Ducks and reexamine the potential role of adaptive evolution of hemoglobin genes in the establishment of its geographic range using multilocus analyses including data from both neutral and functional loci that for the first time involve dense sampling from the Colombian Andes. Our sampling scheme also allowed us to address the question of whether the Colombian population of Ruddy Ducks may be of hybrid origin as suggested by its wide variation in male facial plumage (Fjeldså, 1986;McCracken & Sorenson, 2005;Siegfried, 1976). Our extended F I G U R E 1 Current geographic distribution of Ruddy Ducks and locations where samples were obtained for our phylogeographic analyses sampling suggests that there is likely more complexity to biogeographic and evolutionary scenarios proposed by earlier studies.

| Sampling
To complete the geographic sampling of Ruddy Ducks employed in previous analyses, which only considered a handful of specimens from the Cordillera Central of the Colombian Andes (McCracken & Sorenson, 2005), we obtained 26 additional samples from the Cordillera Oriental of Colombia. These included (1) tissue samples from six specimens from Lake Fúquene (2,580 m), Cundinamarca, obtained from the Museo de Historia Natural at Universidad de los Andes, and (2) 20 blood samples we collected from Ruddy Ducks captured and released at two adjacent artificial lakes (i.e., abandoned gravel pits) located near the town of Guasca, Cundinamarca (2,630-2,640 m elevation; see Table S1).
To examine the distribution of genetic variation within and among North America, Colombia, and the Southern Andes, we are genetically distinct clusters, we used the multilocus clustering algorithm in the program STRUCTURE v 2.2.3 (Pritchard, Stephens, & Donnelly, 2000). We based this analysis on our seven-locus data set, where each locus was coded as a haplotype, and included all the individuals for which no more than two of the loci were missing (n = 85 individuals). We used the admixture model (α = 1) with correlated allele frequencies and tested models for varying numbers of populations (K = 2-4). Each analysis was run for one million generations (burn-in 100,000 generations) and replicated 10 times.
To choose the best value of K, we used the Evanno, Regnaut, and Goudet (2005) method. The outputs of the repeated runs at each K were summarized using the greedy algorithm in CLUMPP v1.1.2 (Jakobsson & Rosenberg, 2007).

| Coalescent analyses of migration
To estimate gene flow (M = m/u) between populations, we used a population genetic approach based on isolation-with-migration coalescent models implemented in the program IMa2 (Hey, 2010). Using a two-population model (i.e., multiple pairwise comparisons), this analysis allowed us to examine the direction of gene exchange based on estimates of historical gene flow between the three populations (Hey, 2005(Hey, , 2010Hey & Nielsen, 2004;Nielsen & Wakeley, 2001).
Because IMa2 assumes no intralocus recombination, we tested for evidence of recombination in our data using the four-gamete test (Hudson & Kaplan, 1985) in DnaSP v.5.10 (Librado & Rozas, 2009). GRIN1 and β A showed evidence of recombination; for the former, we excluded five individuals following Muñoz-Fuentes et al. (2013) and for the latter, we removed the recombinant blocks by eye. Because mtDNA is maternally inherited, we used an inheritance scalar of 0.25 and used the HKY mutation model for this locus (Hasegawa, Kishino, & Yano, 1985). For biparentally inherited nuclear genes, we used an inheritance scalar of 1.0 and implemented the infinite alleles model (Kimura, 1969). Because mtDNA and β A may be under selection associated with the occupation of high-elevation environments, we treated these two genes and the other four nuclear loci separately, and also ran the analysis with all genes combined. IMa2 was initially run with wide priors to explore the sensitivity of the parameters to varying upper bounds. These analyses were then repeated with uniform priors including the entire posterior distribution of each parameter observed in preliminary runs. For each run, we reported values of gene flow corresponding to the 90% highest posterior density (HPD) estimated by IMa2. To ensure that all the parameters converged, autocorrelation and effective sample sizes (ESS) were monitored, and runs were continued until the ESS for each parameter was at least 100.

| RE SULTS
Genetic structure between populations of Ruddy Ducks varied considerably across the seven loci, with Φ st values ranging from 0.12 to 0.61 (Table 1). Despite this variation, genetic differentiation was significant for all population comparisons and loci, except between the Southern Andes and Colombia for the FGB intron.
We found a total of 25 mtDNA haplotypes; 15 of these were present only in North America, five occurred only in the Southern Andes, one was shared between North America and the Southern Andes, one was shared between Colombia and the Southern Andes, and three were shared between North America and Colombia ( Figure 2). Of the four haplotypes found in Colombia, three (including the most common haplotype in North America) were shared with North America, and the other one was the most common haplotype of the Southern Andes.
Varying patterns of genetic variation were observed in nuclear introns ( Figure 3). The most common FGB allele was found in the Southern Andes and also occurred at high frequency in Colombia; six different FGB alleles were found in North America, and two of these were found in Colombia as well. For GRIN1, a private allele and a shared allele were found in Colombia, and the most common allele was shared among the three populations. Of the remaining eight GRIN1 alleles, one was private to the Southern Andes, five were private to North America, one was shared between North America and the Southern Andes, and one was shared between North America and Colombia. Two ODC1 alleles were shared between Colombia and the Southern Andes, and between North America and Colombia; three alleles were private to North America, and one was private to Colombia. The most common PCK1 allele was found in all three populations; one additional allele was private to the Southern Andes, and two were private to North America.
We found eight α A globin alleles, three of which were private to North America, two private to the Southern Andes, one was shared between Colombia and North America, and the two most common alleles were shared by the three populations. Allelic richness standardized to the smallest sample size was greater in North America than in Colombia for all but one locus (ODC1), and more allelic diversity was found in Colombia than in the southern Andes for three of five loci (FGB, ODC1, and HBB). However, mtDNA haplotype diversity was greater in the Southern Andes than in Colombia (Table 2).   were positive (Figure 6g). In turn, gene flow from Colombia to the Southern Andes appeared to be zero, but gene flow into Colombia from the Southern Andes appeared to have occurred (m/u = 1.94; Figure 6h).

| Biogeography of Ruddy Duck populations
A previous phylogeographic analysis was consistent with the hypothesis that Ruddy Ducks experienced a stepwise colonization of the  to the Southern Andes was clearly zero ( Figure S1). In contrast to the scenario proposed by Muñoz-Fuentes et al. (2013), these results may be consistent with a south-to-north colonization of highland areas from the southern South American lowlands, a pattern seen in different waterfowl species (Fjeldså, 1985  sequences between lowland and highland Ruddy Ducks . However, such analyses did not rule out the possibility that genetic variants contributed to variation in other structural or functional properties of the hemoglobin protein . Our results suggest this possibility is unlikely given that we found that not all Colombian individuals have the "high-elevation" genotype identified by the earlier study: only 20% of the individuals from Colombia were homozygous for the allele previously thought to be restricted to high elevations, 52% were heterozygous, and 28% had the lowland genotype found in North America. This result, together with functional analyses  suggest that a Thr/Ser-69 substitution arising de novo did not facilitate the colonization of high-elevation environments by Ruddy Ducks, implying that the previously proposed evolutionary scenario requires reevaluation. Other studies have revealed that patterns of elevational differentiation in hemoglobin sequences consistent with spatially variable selection need not imply adaptive evolution because amino acid replacements may have no effect on respiratory properties of hemoglobin isoforms (Cheviron et al., 2014). However, we note that recent analyses indicate that Ruddy Ducks throughout their range, including North America, possess hemoglobin with relatively high O 2 affinity, similar to other Andean species .

| Hybrid origin of Colombian Ruddy Ducks?
Our results are consistent with the hypothesis suggested by several authors (Fjeldså, 1986;McCracken & Sorenson, 2005;Siegfried, 1976

| CON CLUS IONS
Our work echoes recent suggestions that examining patterns of genetic and phenotypic (i.e., physiological, plumage) variation with dense sampling in populations from northern South America is critical to gain a comprehensive understanding of the biogeographic history of species with broad distribution ranges in the Americas (Avendaño, Arbeláez-Cortés, & Cadena, 2017;Pérez-Emán et al., 2017). Specifically, relative to previous analysis, our addition of samples from the Colombian Andes to phylogeographic analyses revealed (1) that the direction of colonization of Ruddy Ducks appears more difficult to determine-and may even be opposite-than that which was previously inferred, and (2) that if adjustment to high-elevation conditions played an important role in colonization of new areas by Ruddy Ducks, it likely occurred through different adaptive mechanisms than those previously considered or via phenotypic plasticity in physiological parameters that remain to be studied. In any event, our study also exemplifies the breadth of inferences about species history one can potentially make by combining data from nuclear and functionally important loci in phylogeography .

ACK N OWLED G M ENTS
We thank the people whose valuable field assistance helped us to gather data in Colombia, Perú, and in the United States: N.
Wilner. We thank Agregados del Norte for their support during field visits in Colombia. We thank J. Fjeldså for allowing us to use and modify his illustrations. We also thank Autoridad Nacional improved thanks to comments from two anonymous reviewers.

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R CO NTR I B UTI O N
MLJ, CDC, and KMC concieved the ideas and designed the methodology. MLJ and KMC conducted fieldwork and collected data. MLJ analyzed the data. MLJ and CDC led the writing of the manuscript.
All authors contributed critically to the drafts and gave final approval for publication.