Introduction

The role of lowland habitats in the accumulation of tropical diversity has been the subject of considerable research (Mittelbach et al. 2007, Rull 2011). A great deal of scrutiny has been directed upon Amazonia in particular and the patterns of endemism found therein (e.g. Tuomisto et al. 1995, Wesselingh et al. 2009, Hoorn et al. 2010, Cheng et al. 2013). However, studies from complementary lowland biomes are necessary to develop a complete picture of habitat change in the low-elevation tropics as a whole. The evolutionary history of lowland tropical savannas and seasonally dry forests in the Americas is essential for an integrated understanding of biome dynamics across space and time. We use both microsatellite markers and mitochondrial DNA sequence to document phylogeographic and population-level genetic structure in the burnished-buff tanager (Tangara cayana Linnaeus), a species that is ubiquitous and widespread throughout South American savanna and seasonally dry forests. We relate these results to inferred T. cayana distributions during the mid-Holocene and Last Glacial Maximum (LGM) using environmental niche modeling.

The biogeography of lowland open biomes in South America is complex, and our understanding of the evolutionary dynamics for this suite of ecosystems is far from complete (Werneck et al. 2011). The open woodlands, savannas, and tropical dry forests of South America are patchily distributed along the periphery of Amazonia (Fig. 1), and natural histories can vary from one area to another (Linares-Palomino et al. 2011). These habitats include seasonally dry forest nuclei (e.g. Caatinga), tropical savannas and woodlands (e.g. Cerrado, llanos Pampas del Heath, Pantanal), and several inter-montane valleys along the Andes in western South America (e.g. Huallaga and Mayo River valleys in Peru). The complexity of the various open lowland biomes in the Tropical Americas reflects their disparate evolutionary histories. Some authors have posited that the circum-Amazonian distribution is the result of a “dry forest refugia” scenario (Prado & Gibbs 1993, Pennington et al. 2000), in which formerly extensive tracts of dry forest have contracted toward the present. An alternate scenario is that these habitats have always been fragmentary and isolated, with floristic and faunistic similarities resulting from dispersal (Mayle 2004, 2006).

We study the intraspecific genetic structure of a species with a typical “circum-Amazonian” distribution, T. cayana. The patchy nature of T. cayana's distribution is ideal for examining the competing hypotheses of vicariant refugia versus long-distance dispersal. Shared haplotypes among populations may indicate that gene flow occurred during the recent past. Conversely, populations which have coalesced may be consistent with relatively rare colonization events via long-distance dispersal with no subsequent gene flow.

The implications of these competing hypotheses for the historical biogeography of seasonally dry tropical forests are profound. Evidence in support of vicariant refugia implies large-scale expansion of suitable T. cayana habitat over much of South America. Alternately, support for long-distance dispersal may indicate that tropical dry forest centers are quite stable, undergoing relatively modest episodes of expansion and contraction.

T. cayana is a member of the largest genus of passerine birds in South America. Species of the genus occur in nearly all forested biomes of Central America and northern South America from the Atlantic and Amazonian forests in the lowlands up to about 3500 m in the cloud forests of the Andes. The systematics and phylogeny of the Tangara have been studied in detail (Isler & Isler 1999, Burns & Naoki 2004, Sedano & Burns 2010, Burns et al. 2014). T. cayana is a member of a clade of seven allospecies spread throughout continental South America, Central America, and the Lesser Antilles. The monophyly of the clade has been strongly supported in molecular studies (Sedano & Burns 2010, Burns et al. 2014). Of the seven, T. cayana has by far the broadest geographic distribution despite being one of the younger members of the clade (Burns & Naoki 2004). Other members of this clade include rufous-winged tanager, Tangara lavinia (Central America and the Chocó), black-backed tanager, Tangara peruviana (southern Atlantic Forest), Lesser Antillean tanager, Tangara cucullata (Lesser Antilles), chestnut-backed tanager, Tangara preciosa (southern Atlantic Forest), scrub tanager Tangara vitriolina (intermontane valleys of Colombia) and the greencapped tanager, Tangara meyerdeschauenseei (range restricted in central Andes of Peru). This species group is ecologically distinct from other Tangara species with a preference for open, mostly drier habitats (Isler & Isler 1999). The majority of the clade tends to have relatively restricted distributions, with T. preciosa, T. peruviana, T. cucullata, and T. meyerdeschauenseei all being examples of tanagers with particularly narrow geographical ranges. Within this group, T. cayana is the notable exception, occurring in the tropical Savannas that are widespread both to the north (“llanos”) and south (“cerrado”) of the Amazon basin. In addition, T. cayana is patchily distributed to the west of the Amazon in drier patches at the base of the Andes such as the Mayo River valley and the Pampas del Heath as well as in isolated patches of savanna habitat within Amazonia. Together, this range forms a classic “circum-Amazonian” distribution (Remsen et al. 1991, Bates 1997).

The origin of the circum-Amazonian distribution has been the subject of some speculation. Although longdistance dispersal cannot be ruled out, the patchiness of the present T. cayana distribution to the west of the Amazon suggests that its range may have been more extensive in the past (Mayle 2006). Populations in areas that are now isolated by distances of thousands of kilometers were likely established at a time when suitable habitat for T. cayana was more contiguous.

T. cayana has seven described subspecies and considerable phenotypic variation. However, differences are frequently clinal despite T. cayana's disjunct distribution, with the most divergent phenotypes found in the northern and southern extreme of the species' range (Isler & Isler 1999).

Material and Methods

To model the current and historical distributions of T. cayana, we used the program MaxEnt (Phillips & Dudík 2008), which produces a map of continuous habitat suitability scores. To model distributions, we gathered locality data from museums, accessed via the GBIF portal. Duplicate occurrence records were trimmed using the program ENMTools (Warren et al. 2010), as were any occurrences that fell over 100 km outside of the known documented distribution of T. cayana. In addition, we used the ArcGIS Toolkit SDMtoolbox (v1.1c, Brown 2014) to rarify spatially clustered observations within a 10 km radius which have the potential to bias model development. Coordinates that fell within the same pixel in the Worldclim raster were likewise trimmed for a total of 434 unique localities. Localities of tissues used in subsequent genetic studies were not included, because many of them lacked the precise georeferences characteristic of the MaxEnt dataset. To represent climatic conditions, we used the Worldclim climatic layers (Hijmans et al. 2005) at 2.5 arc-minute resolution and cropped to 15° N, 60° S, 90° W and 30° W to encompass all of South America. To avoid overfitting of the environmental niche models, highly correlated Worldclim layers were eliminated. In total, six layers were selected based on uniqueness and relative importance in model development. The aforementioned six Worldclim layers are detailed in Table 1. Environmental niche models were projected onto two alternate climate models at the LGM ∼21kybp, and at the mid-Holocene, ∼6kybp. These alternate paleoclimate reconstructions are based upon the general circulation models MIROC-ESM and CCSM4 (Otto-Bliesner et al. 2006). In estimating the environmental niche model, twenty replicates were developed and their performance was assessed via cross-validation. We used the overprediction correction function in SDMtoolbox (v1.1c, Brown 2014) to crop the layers to a polygon corresponding to the geographic limits of T. cayana occurrences used in the model development, with a 100 km buffer zone.

Fig. 1.

Geographic localities of individuals included in this study, along with the distribution of the Seasonally Dry Tropical Forest (SDTF) and Savanna biomes (modified from Werneck et al. 2011) that largely define the distribution of Tangara cayana. Colors of the sampling localities correspond to the approximate microsatellite group assignment of those populations as determined by STRUCTURE analysis. To the right is the STRUCTURE analysis of T. cayana microsatellites indicating approximate individual genotype assignment to three microsatellite groups. Inividuals belonging to the red microsatellite group were recovered in two different localities in Pernambuco Brazil along with admixed individuals. The same microsatellite group is also found admixed in populations north of the River Amazon.

Table 1.

BIOCLIM layers used in development of the ENM for Tangara cayana, and their respective definitions.

Genetic data

Extraction of whole DNA from tissue samples was performed using a Qiagen DNeasy kit (Qiagen U.S.A., Valencia, CA) using the standard protocol as recommended by the manufacturer for animal tissue. For analysis of phylogeographic and demographic patterns of intraspecific divergence, we gathered data from both nuclear (microsatellite) and mitochondrial sources. The mitochondrial data gathered consisted of cytochrome b and the cytochrome oxidase three (CO3) genes sequence (maximum coverage of 1020 and 715, respectively). Amplification via PCR was performed using the following protocol for both genes: after an initial denaturation step at 96 °C for five minutes, 40 cycles were performed. Each cycle consisted of a denaturation step at 96 °C for 30 seconds, an annealing temperature of 50 °C for 30 seconds, an extension temperature of 72 °C for one minute. At the end of the 35 cycles, the PCR product was held at 4 °C indefinitely. Amplification protocol was nearly identical for microsatellite loci, with the only exception that there were 50 cycles instead of 40, and annealing temperature was gradually ramped up from 50 °C to 65 °C over the course of the 50 cycles. PCR product was cycle sequenced using a BigDye Terminator 3.1 cycle sequencing kit (Catalog #4337455 Life Technologies, Grand Island, NY), cycle sequenced for 35 repetitions at an annealing temperature of 50 °C. Samples were sequenced on a 3730xl DNA Analyzer (Applied Biosystems, Carlsbad, CA). Sequences were corrected and edited manually in Sequencher (version 5.1, Gene Codes Corporation, Ann Arbor, MI), and Mega v5.1 (Tamura et al. 2011).

Table 2.

Sampling localities for Tangara cayana tissues analyzed in this study. Abbreviations for tissue sources are as follows: FMNH = Field Museum of Natural History, MPEG = Museu Paraense Emílio Goeldi, LSUMNS = Louisiana State University Museum of Natural Science, ANSP = Philadelphia Academy of Natural Science.

The microsatellite data set consists of 10 loci, the primers for which were developed by Corrêa et al. (2010) for the thraupid species Neothraupis fasciata, and which were shown to successfully amplify across a wide suite of thraupid taxa. All of these loci were polymorphic within and/or between individuals. The use of microsatellite primers designed for an outgroup avoids any potential for development bias toward a particular Tangara species. This same lack of development bias may also result in the use of microsatellite markers that are less variable than those that could be developed for a particular species. However, given the time scale that separates many populations, highly variable markers would be of limited utility and potentially confounding. Microsatellites were amplified via PCR using fluorescent-tagged nucleic acids (Schuelke 2000) using the protocol described above. Fragment length was measured using a 3730xl DNA Analyzer (Applied Biosystems, Carlsbad, CA), and results were processed using Genepop software (version 1.2, Raymond & Rousset 1995).

Microsatellite analysis

We performed a STRUCTURE analysis on the microsatellite data set (version 2.3, Pritchard et al. 2000, Hubisz et al. 2009). This algorithm has the advantage that population membership is not specified a priori. Rather, probability of group membership is assigned a posteriori. The first 100000 iterations of each MCMC chain were discarded as burn-in, and an additional 1000000 steps were logged thereafter. For each potential number of microsatellite groups (k) from one through ten, we performed ten replicates. Output from STRUCTURE analyses was processed using STRUCTURE HARVESTER (version 0.6.93, Earl & vonHoldt 2012), and further analyzed using CLUMPP (version 1.1.2, Jakobsson & Rosenberg 2007). We utilized the Evanno method as implemented in STRUCTURE HARVESTER to identify the number of microsatellite groups with the highest likelihood. Graphs of the results of the STRUCTURE analysis were produced using the program Distruct (Rosenberg 2004).

Fig. 2.

Environmental niche model for Tangara cayana projected onto present day bioclimatic variables (left), as well as two alternate models of paleoclimate at the Last Glacial Maximum (CCSM, center and MIROC, right). Niche models are projected onto paleoclimate layers at both Last Glacial Maximum (∼22 kybp, top) and mid-Holocene (∼6 kybp, bottom). Warmer colors indicate a greater degree of habitat suitability. Black outline on the present-day ENM is a published range map for this species (Birdlife International & Nature serve 2012). Note that this range map omits populations to the west of the Amazon in Peru and Bolivia which are documented elsewhere (e.g. Isler & Isler 1999). Pink dots represent T. cayana localities used in model development.

Sequence analyses

For analysis, sequences from the two mitochondrial genes were concatenated, since they are nonrecombinant and evolve as a unit. We estimated ultrametric phylogenies in BEAST using the HKY substitution model a relaxed clock (Drummond et al. 2006), with 10 million iterations of the MCMC chain and the first 10 % of trees discarded as burnin. Results of the BEAST analyses were processed in TreeAnnotator (version 1.7, Drummond et al. 2012) to obtain a consensus tree and posterior probabilities for the nodes. To root the tree, we used data from Tangara chilensis, as sequences for the CO3 gene were not available from closer relatives to T. cayana. To estimate divergence times we used the ultrametric phylogeny and assumed a substitution rate of 2 % per million years (Lovette 2004a, Päckert et al. 2007). We used the program Tracer (version 1.7, Rambaut et al. 2014) for parameter analysis to ensure an effective sample size of > 200 for all parameters.

Fig. 3.

Intraspecific phylogeny for Tangara cayana. Asterisks (*) indicate nodes with posterior probability > 0.98. Clade colors refer to geographic origin as follows: blue = southeast Brazil, green = Bolivia, black = Guyana and red = northeast Brazil. Scale bar indicates 0.5 % uncorrected sequence divergence.

We constructed a Median-Joining haplotype network (Bandelt et al. 1999) using the program PopART (Leigh 2015). The same program was also used to produce the graphical representation of the haplotype network.

Results

Environmental niche modelling and paleoclimate projection

The Area Under the Curve (AUC, a measure of model performance) of the replicate MaxEnt runs was 0.868 with a standard deviation of 0.018. The species distribution model for the present is consistent with published accounts of the species range (Fig. 2). Both paleoclimate models indicate more suitable habitats for T. cayana at mid-Holocene and LGM, and agree broadly on several points. Both MIROC and CCSM show increased suitable habitats along the base of the Andes particularly in southwest Amazonia. This habitat suitability was highest in both models at LGM and in both cases decreased toward the present, though CCSM model shows a pattern much closer to presentday values by the mid-Holocene. Both models also indicate more suitable T. cayana habitat during LGM across much of the Guyanas, Rondonia, and northeast Amazonia. Both paleoclimate models indicate large, contiguous areas of suitable habitat in western Amazonia toward the LGM. While both models show good habitat suitability throughout the Atlantic forest region, there are also large swathes of adjacent habitat in southeastern Amazonia that is unsuitable for T. cayana from the present through the LGM. The distribution of suitable habitat in these areas varies between the two paleoclimate reconstructions in these areas. Both paleoclimate models also indicate varying degrees of low habitat suitability in southeastern Amazonia during the LGM.

Fig. 4.

Median Joining haplotype network for Tangara cayana (Bandelt et al. 1999). Size of circles indicates the number of individuals found with that haplotype. Small black circles represent hypothetical haplotypes. Hash marks are equivalent to one base mutation. Coors correspond to sampling location and geographic group as in Fig. 3. Numbers refer to accession IDs listed in Table 2.

Discussion

T. cayana is a comparatively young species for the genus, having diverged from the Lesser Antillean tanager (T. cucullata) within the last 800000 years according to a previously published phylogeny (Sedano & Burns 2010). Levels of uncorrected sequence divergence in the intraspecific T. cayana phylogeny in this study are as high as 1 %, or 500000 years assuming a rate of 2 % sequence divergence per million years (Lovette 2004a, Päckert et al. 2007). This is in concordance with previously published genus-level phylogenies (Burns & Naoki 2004, Sedano & Burns 2010, Burns et al. 2014).

The most closely related taxa to T. cayana are T. cucullata and T. vitriolina, which are found in the Lesser Antilles and dry intermontane valleys of Colombia, respectively (Burns et al. 2014). However, the basal divergence in the intraspecific T. cayana phylogeny presented here subtends the populations from southeastern Brazil including the states of Piauí, Pernambuco, and Bahia (Fig. 3), though it should be noted that this node has very weak support in the Bayesian phylogeny. The next most divergent clade within T. cayana corresponds to the populations in Bolivia. The remaining populations, from north of the River Amazon, fall into two groups. Individuals from Guyana are a strongly supported clade which is basal to the rest of the northern populations. The remaining populations east of the Guyanas in northeastern Brazil show little discernable structure. The phylogeographic structure recovered in T. cayana suggests an origin in the south. Subsequently there was expansion from Bolivia probably along the base of the Andes to reach Guyana and finally northeastern Brazil. The potential relationship between Guyanan and northern Amazonian populations is observed on the Bayesian phylogeny, but not in the haplotype network analysis. The recovery of a highly divergent “Bolivian” haplotype in the Guyanan population is consistent with a historical connection across western Amazonia, as predicted by ENM projections onto both paleoclimate models. Environmental niche models also support the hypothesis of a colonization route through western Amazonia along the base of the Andes. Both CCSM and MIROC paleoclimate models show greater degrees of habitat suitability and connectivity during the LGM and decreasing through the mid-Holocene toward the present, particularly in western Amazonia (Fig. 2).

Mitochondrial sequence and microsatellite data show generally congruent patterns, although there is evidence for some additional diversity in the microsatellite data. The mitochondrial phylogeny groups the populations from Taquaritinga in Pernambuco state, Brazil, with the basal clade including populations from Bahia and Piauí (corresponding to the blue microsatellite group in Fig. 3). Other populations in Pernambuco state are grouped with a clade that includes populations to the north of the Amazon. In the microsatellite data, these Pernambuco populations are assigned to a different microsatellite group (red in Fig. 1) which predominates in Pernambuco populations, but which is also found in T. cayana north of the Amazon. Thus, both microsatellite and mtDNA sequence data support a relationship between eastern Pernambuco populations and those north of the River Amazon. The microsatellite data further suggest recent or ongoing admixture between northern populations and those in eastern Pernambuco (Fig. 2). While the mitochondrial phylogeny has a strongly supported node that separates Bolivian T. cayana from those north of the Amazon, the microsatellite data largely assigns these populations to the same microsatellite group (green in Fig. 2). The Bolivian populations show no evidence of admixture; however, populations east of the Guyanas in northern Brazil show increasing evidence of admixture with eastern Pernambuco populations. Relative to the northern Brazilian birds, the Guyanan samples have a lower degree of admixture with the Caatinga microsatellite group, and their microsatellite profile is comparable to Bolivian T. cayana (Fig. 1). An affinity between the Guyanan samples and those from Bolivia is further evidenced by the mitochondrial data, in which a Bolivian haplotype clusters with that population (see Fig. 3). The co-occurrence of two otherwise geographically segregated haplotypes in Guyana suggests it as an area that has had some degree of gene flow from different biogeographic regions during the last several hundred thousand years. However, the Bolivian haplotype recovered from the Guyanas is consistent with introgression from the Bolivian population into the Guyanas and/or possible retention of ancestral polymorphism.

Microsatellite data indicate two distinct groups centered around the Caatinga and other seasonally dry forests in southeastern Brazil (Fig. 2). The geographic ranges of these microsatellite groups closely coincide with areas of endemism that have been documented across a variety of vertebrate taxa (Carnaval et al. 2009). This includes anurans (Carnaval & Bates 2007, Carnaval & Moritz 2008) and birds (Cracraft 1985, defined therein as the “Parana Center”). While there has historically been a paucity of documented endemism among mammals (Mares et al. 1985), recent work has documented considerable differentiation in the tigrina (Leopardus tigrinus) between Caatinga populations and those further to the south (Trigo et al. 2013).

Many avian species and species groups have a “circum-Amazonian” distribution (sensu Remsen et al. 1991). This distribution is characterized by populations distributed along the base of the Andes, as well as north and south of the Amazon. The distributions of dispersal-prone volant species suggest at least two different evolutionary scenarios for seasonally dry forest species. Often, these species lack a clear eastern connection on the Atlantic coast, particularly around the Amazon delta. However, there are many species which show a conflicting distribution, with largely continuous ranges to the east and a lack of connectivity along the base of the Andes (e.g. da Silva & Bates 2002). Previous studies examining the evolution of the circum-Amazonian distribution have shown ambiguous and sometimes conflicting patterns. For example, Lovette (2004b) finds a strongly supported relationship between Guyanan and Atlantic forest populations in the Phaeothlypis wood-warbler complex, while a similar study of the Synallaxis ruficapilla species group showed evidence for an affinity between Caatinga populations and those in the N Andean foothills (Batalha-Filho et al. 2013). The possibility of a connection between Guyanan and Bolivian populations via the “dry forest arc” along the base of the Andes is a pattern that has also been observed in other taxa (e.g. rattlesnakes, Wüster et al. 2005).

Studies of understory taxa with circum-Amazonian distributions have recovered considerable levels of divergence on the order of 1–2 million years (Lovette 2004a, Batalha-Filho et al. 2013), which far exceeds the divergence in T. cayana documented here. While many species with circum-Amazonian distributions also show considerable morphological variation (Bates 1997), this variation may not necessarily coincide with genetic breaks between populations (Irwin 2002, Lovette 2004b). Previous studies have established that ecology is directly related to dispersal ability (e.g. Terborgh 1974, Faaborg 1979, Lewis et al. 2004, Moore et al. 2008), and that this in turn is a major predictor for the amount of differentiation between populations or within a taxon (Burney 2009, Smith et al. 2014).

Foraging stratum has also been shown to have an important effect on lineage diversification (Smith et al. 2014). Tangara in general are canopy species, and canopy dwellers tend to have lower levels of intraspecific differentiation relative to understory species (Burney & Brumfield 2009). T. cayana inhabits fairly open habitats with patchy tree cover, and thus may be even more prone than many other Tangara species to disperse over fairly long distances. Bonaccorso et al. (2006) presented evidence from environmental niche modeling which suggests that species inhabiting drier habitats around the periphery of Amazonia may have had varied responses to climate fluctuations. As such, the timing and route of T. cayana's spread across its circum-Amazonian distribution may not be applicable to many other circum-Amazonian taxa. These disparities among particular evolutionary histories suggests that current circum-Amazonian distributions may have come to be established several different ways.

Conclusions

The evolution and diversification of T. cayana across the breadth of its geographic range informs our understanding of seasonally dry forest and forested savannas across time and space. The mtDNA data for T. cayana illustrate a pattern of dispersal and diversification involving a likely origin in southeastern South America, with subsequent colonization of the circum-Amazonian dry forests in a clockwise fashion; beginning in the south, expanding west and proceeding northward along the base of the Andes toward the Guyanas, finally crossing into northeast Brazil. The presence of a southern (Bolivian) haplotype in the Guyanan population may be indicative of historical connectivity across the dry forest arc distributed along the base of the Andes. The position of the Guyanan samples at the base of the Bolivian clade indicates that this haplotype may reflect retention of ancestral diversity and is not necessarily a result of recent introgression or gene flow. The Guyanan populations of T. cayana are crucial to understanding the evolutionary history of the species and warrant further study. There is also evidence in the microsatellite data for a relatively recent introgression of northern Brazilian populations into eastern Pernambuco to the south of the Amazon. Data from mtDNA sequence and microsatellites are in broad agreement, both supporting discrete groups found in southeast Brazil and an affinity between Bolivian populations and those to the north of the Amazon. Microsatellite data indicate additional diversity in the Caatinga populations in Pernambuco. This diversity is congruent with an area of endemism based on differentiation within a wide variety of vertebrate taxa. Microsatellites do not show a clear distinction between Bolivian populations and those to the north of Amazonia. However, populations to the north are more highly admixed with the “Caatinga” microsatellite group moving further east.

While the available data suggest that this species successfully colonized patches of suitable habitat via a western route, better sampling of the highly isolated T. cayana populations in the Andean dry forest arc is needed. In particular, the inclusion of samples from the two isolated populations in Peru (Pampas del Heath and the Mayo Valley) would help to clarify the timing and precise route by which the Andean dry forest arc was colonized by T. cayana.