Late Pleistocene climatic changes promoted demographic expansion and population reconnection of a Neotropical savanna-adapted bird, Neothraupis fasciata (Aves: Thraupidae)

We performed phylogeographic and genetic structure analyses of Neothraupis fasciata joined with species distribution modelling to evaluate whether: (1) the distribution of genetic variability shows a pattern expected by the isolation-by-distance model; (2) the influence of the Pleistocene climate changes on species distribution; and (3) climate/climatic stability (hypothesis of climatic stability) as a predictor of population genetic diversity. Based on two molecular datasets (ND2 and FIB-5), the isolation-by-distance hypothesis was not supported. The mitochondrial haplotype network indicated the existence of historically isolated populations at the southern range of the species distribution, and recent population expansion was identified by both neutrality tests and extended Bayesian skyline plot analysis. Thus, the climatic changes during the Pleistocene might have promoted the reconnection of the partially isolated southern populations, which may have persisted in the plateaus during the cycles of savanna contractions. Subsequently, this species (re)colonized northern areas of the species present distribution, following the continuous vegetation on the São Francisco and Central plateaus about 60 kyr, and also reached the Amazonian savannas likely via the central corridor. Thus, our results indicated that the intrinsic relationship between the relief heterogeneity (plateaus and depressions) and the climatic fluctuations, mainly in the Pleistocene, promoted population reconnection and demographic expansion of N. fasciata.


Introduction
Tertiary geological events and the Quaternary climatic changes have been considered the key events related to the biodiversity diversification (e.g. [1][2][3][4][5] climatic cycles had great importance in South American forest and nonforest vegetation (e.g. Savannas) range dynamics [1,3,[6][7][8][9][10][11][12]. Briefly, it is postulated that during moist and warm periods (interglacial periods), forests expanded their ranges, while during the cold and dry periods (glacial periods) the savannas reached their maximum extension (see [9]). Events of both range contraction and expansion during glacial periods have been hypothesized for the "morphoclimatic domain of the Cerrado" [13,14], hereafter referred to as the Cerrado. For instance, during the glacial period the current range of this savanna was replaced by Araucaria forests at its southern edge and by xeric forests at its northern edge, while it was expanded throughout eastern and central Amazonia [14]. According to Werneck et al. [3], during the Last Interglacial (LIG, c. 120,000 years ago or 120 kyr) the Cerrado expanded, while during the Last Glacial Maximum (LGM, c. 21 kyr) the Cerrado reached its smallest range. It is also postulated that during the Pleistocene the forest and nonforest vegetation range dynamics may have led to the establishment of biogeographic corridors currently connecting the isolated savanna blocks in northern Amazonia Forest (Llanos and Amazonian savannas) and Cerrado [1,9,[15][16][17][18][19][20]. The dynamic of range shifts of the Cerrado is complex and there is still no clear pattern related to the diversification of the Cerrado biodiversity that has emerged driven by Pleistocene climatic fluctuations. Considering that the Cerrado is one of the biodiversity hotspots [21][22][23], understanding the processes that led to the diversification of its biodiversity is also fundamental for its conservation. The demographic history of the South American species occurring in open vegetation formations was probably influenced by the climatic fluctuations, with expected events of population expansion (see leading-edge model [24,25]) following the expansion of suitable habitats and population bottlenecks or extinctions (see rear-edge model [24,25]), in response to the retraction of habitat suitability [3,11]. Cycles of population expansion and/or contraction may promote secondary contact among allopatric populations or different episodes of colonization, resulting in different genetic signals. For events of secondary contact, an increase in genetic diversity in the population that received more migrants is expected, since migrants are likely to carry new alleles to the receiver population (see [26]). On the other hand, colonization events from a core population lead to a gradual reduction in genetic diversity towards newly colonized areas due to genetic drift operating through repeated founder events or bottlenecks [24][25][26][27]. Additionally, when the equilibrium between migration and genetic drift is reached, patterns of genetic diversity may be also explained by the isolation-by-distance model, in which geographically close individuals tend to be genetically more similar than individuals that are further apart [28,29]. Finally, areas that remained climatically stable through the Quaternary climatic changes are expected to exhibit higher genetic diversity than unstable ones, since unstable areas must have experienced more extinction and recolonization events than stable areas [30,31].
Herein, we performed phylogeographic and population genetic analysis of Neothraupis fasciata (Lichtenstein, 1823) to evaluate the response of this savanna-adapted bird to environmental changes during the late Quaternary (last 120 kyr). Commonly known as a Whitebanded Tanager, it is a passerine species from the Thraupidae family, and is the only representative of the genus. This tanager is widely distributed and very common in the undisturbed Cerrado, occurring mainly in cerrado sensu stricto and grassy cerrado [32][33][34][35]. This species is resident and territorial [36,37], has a low dispersal capacity, with a mean dispersal distance of 200 meters per breeding season [37,38]. A previous study using a set of microsatellite markers found moderate genetic structure in N. fasciata and the distribution of the genetic variability was not explained by the isolation-by-distance hypothesis, habitat heterogeneity or by the core-periphery effect [39]. Similarly, the intense biome fragmentation by anthropic activity has not contributed significantly to the current pattern of genetic structure found in this species [39]. These findings suggest the influence of historical processes on the distribution of genetic diversity of this tanager.
Based on molecular data sets and ecological niche modeling, we evaluated: (1) the isolationby-distance hypothesis, with an increase in genetic differentiation among locations being expected due to an increase in geographic distances; (2) the influence of late Quaternary climate changes on the species distribution and demographic history, with events of population expansion being expected in response to the increase in suitable areas and a reduction in genetic diversity following population expansion; and (3) climate/climatic stability (hypothesis of climatic stability) as a predictor of population genetic diversity, with greater genetic diversity in stable areas being expected than in unstable ones.

Species distribution modeling
The species distribution models were estimated compiling occurrence records of N. fasciata from fieldwork of the Laboratório de Genética e Biodiversidade da Universidade de Brasília, scientific articles, museums, and using the two online databases (Global Biodiversity Information Facility: http://www.gbif.org/; SpeciesLink: http://splink.cria.org.br/). We filtered the 176 unique occurrence points found by excluding redundant records in each pixel on the scale of the bioclimatic variables, resulting in 158 points. A set of 19 bioclimatic variables and elevation for the present were downloaded from WorldClim website in a 2.5 arc-min resolution (http:// www.worldclim.org/). The highly correlated variables (|r|�0.75) were eliminated, and variables with |r|<0.75 were retained using findCorrelation function available on caret (Classification And Regression Training) package [40] in R 3.4.2 [41]. The correlation between the retained bioclimatic variables (Bio2, Bio3, Bio10, Bio13, Bio15, Bio18, and Bio19) ranged from -0.579 (Bio10~Bio18) to 0.531 (Bio3~Bio 15).
The maximum entropy approach implemented in MAXENT 3.4.1 [42] was used to build distribution models [43]. Distribution models were run with 20 replicates subsample, using 20% of the points for test, 10,000 background points, 1,500 maximum iterations and other parameters set as default. Discrimination capacity models were evaluated based on the average values of the area under the curve (AUC), sensitivity, specificity and accuracy. We used a set of seven non-correlated bioclimatic variables to estimate the species distribution under mid-Holocene (6 kyr), LGM, and LIG climate scenarios. Past conditions were based on two general circulation models (CCSM4 and MIROC-ESM) for Holocene and LGM, and the data simulated by Otto-Bliesner et al. [44] for the LIG. All bioclimatic datasets were downloaded in a 2.5 arc-min resolution, except for LIG, which was downloaded at 30 arc-sec resolution and exported to a 2.5 arc-min resolution by calculating the mean value for each variable using the aggregate tool in ARCMAP 10.1 (Esri).
To convert habitat suitability into presence/absence binary maps, we used the equal training sensitivity and specificity threshold. For Holocene and LGM we combined the two binary maps obtained with the CCSM4 and MIROC-ESM general circulation models into a single map, in which only areas with predicted occurrence in both projections were highlighted as potential distribution area. Stable areas over the last 120 kyr were delimited based on binary maps and were defined as areas with potential distribution concordant in all projected periods.
To explore the relief, we used the Global Multi-resolution Terrain Elevation Data 2010 (GMTED 2010 -https://earthexplorer.usgs.gov/), which provides a high level of detail in global topographic data 30 arc-sec, suitable for various regional and continental applications [45]. From the mean elevation data (~1km in the equator), we derived the slope (percent) and minimum curvature (˚.m-1) of Brazil using the System for Automated Geoscientific Analyses (SAGA) software, version 6.2 [46]. We used the elevation, slope and minimum curvature as input grids to run K-mean clustering based on the combined minimum distance [47] and Hill-climbing [48] algorithm, aiming to divide two geomorphological compartments: plateaus and depressions. Considering the habitat preference of N. fasciata, we evaluated the continuity of the plateaus, based on the categories by Riitters et al. [49,50] and the geomorphologic features to delimit the historical population by digital vectorization on a computer screen.

Genetic sampling
Tissue samples were collected from 105 N. fasciata individuals (Fig 1B and 1C) [51], while the main South American savannas are shown in gray [52]. The green arrows indicate the savanna connections (revised by [1]) (a) along the Andes region [20], (b) following the Madeira River [18], (c) directly across central Amazonia [15] and (d) along the Atlantic coast [15].  Table 1). In these areas, individuals were attracted by playback and captured using mist-nets. Blood samples (50-100 μL) were taken from each individual by pricking the brachial vein with a sterile needle and collected with a microcapillary tube, then stored in absolute ethanol. All captured birds were ringed with standard ornithological metal rings supplied by Centro Nacional de Pesquisa e Conservação de Aves Silvestres (CEMAVE) and afterwards released at the capture site. Seven ethanol-preserved muscle from individuals sampled at Brasilândia de Minas (BRMI) and Uruçuí (URUC) were kindly provided by Museu Paraense Emílio Goeldi (Table 1). All work was conducted under the approval of the Ethics Committee of the Universidade de Brasília (UnBDoc n˚75111/2013) and research permits were issued by the Brazilian Environmental Agency Instituto Chico Mendes de Conservação (ICMBio) (SISBIO n˚27682-1).
Total DNA was extracted from blood samples using proteinase K digestion followed by phenol-chloroform extraction [53], while for museum samples, total DNA was extracted using a PureLink Genomic DNA Kit. We amplified the mitochondrial NADH dehydrogenase subunit 2 gene (ND2) by PCR, using primers LMET (J. Groth apud [54]) and H6313 [55] and a fragment of the nuclear beta-fibrinogen gene intron 5 (FIB-5) using the primers FIB5 and FIB6 [56]. The ND2 and FIB introns are markers commonly used in phylogeographic studies of birds and have proven to be informative. PCR reactions were performed in a 10 μl volume containing 30 ng of template DNA, 1X PCR buffer, 0.25 mM of each dNTP, 1 mM of each primer, 1.5 mM MgCl2, and 0.5 unit of Taq polymerase (Invitrogen). Thermocycle conditions were set as: 95˚C for 7 min, then 35 cycles at 95˚C for 1 min / 52˚C for 45 s / 72˚C for 1 min, followed by a final extension at 72˚C for 10 min. We used the Shrimp Alkaline Phosphatase and Exonuclease I protocol for cleanup PCR products according to the manufacturer's instructions. Sequencing reactions were performed with BigDye Terminator v3.1 Cycle Sequencing Kit according to the manufacturer's instructions and PCR products were sequenced in both We edited the sequences in GENEIOUS 6.0.6 (Geneious Co., Wellington, New Zealand) and coded double peaks in both strands of FIB-5 sequence electropherograms using the IUPAC nucleotide code ambiguous positions. In the FIB-5 sequences, we identified alleles with different sizes (14 bp), generating sequences with mixed traces [57] in heterozygous individuals. To solve this sequencing problem, the readable part of the forward sequence trace in which the pair of allelic sequences properly aligned was concatenated with the complement and reverse of the readable part of the reverse sequence trace. The resulting sequence was aligned with the dataset containing homozygous sequences for the indel and the gap between forward and reverse traces was completed with the IUPAC code for any base "N". The alignments of consensus sequences of each dataset were done using the Muscle algorithm implemented in GENEIOUS. For FIB-5 we resolved the gametic phase using the algorithm Phase implemented in DnaSP 5.10 [58] with default configurations.

Molecular analyses
The number of polymorphic sites, number of haplotypes, number of unique haplotypes, haplotype and nucleotide diversities were estimated using DnaSP (Table 1). We estimated genetic diversity indices considering the ND2 and FIB-5 gene fragments separately for each sample site. We constructed haplotype networks for each gene fragment using the median-joining algorithm [59] implemented in POPART 1.7 [60] with all parameters set to default. For FIB-5 network we used the phased dataset inferred in DnaSP.
Genetic differentiation between all sample sites was estimated using F-statistics (F ST ) in ARLEQUIN 3.5.1.3 [61]. Calculations of F ST were done using Tamura-Nei model for ND2 alignment and Jukes-Cantor model for FIB-5 alignment. The best-fit model of evolution of each dataset was selected according to the Bayesian Inference Criterion in JMODELTEST 2.1.10 [62]. We used a total of 50,000 permutations to assess the statistical significance of each pairwise comparison. Genetic groups across the landscape were identified employing a Bayesian approach. The Bayesian Analysis of Population Structure (BAPS) was performed at the population level (spatial clustering of groups), considering each locality as a population and using their respective geographic coordinates in BAPS 6.0 program [63]. To determine the optimal number of genetic groups, BAPS analysis was performed for each marker with the maximum number of groups (K) set to 12, and 20 replicates for each K-value.
To verify whether the genetic differentiation among each pair of individuals sampled in the Cerrado is associated with the geographic distance, we performed a Mantel test. The Euclidean geographic distances were obtained using an on-line distance calculator from the Instituto Brasileiro de Pesquisas Espaciais (http://www.dpi.inpe.br/calcula/). We conducted the test using the corrected genetic distance (PhiST/(1-PhiST)) and the logarithm of the geographic distance in IBDWS 3.23 [64], considering each locus separately and 30,000 randomizations. We also evaluated the influence of geographic gradients on genetic diversity performing a regression analysis of genetic diversity indices (percentage of nucleotide diversity and haplotype diversity) and geographic variables of all sample locations (latitude and longitude, expressed as decimal degrees). All regressions were done using Spearman's rank-order correlation coefficient in PAST 3.11 [65].
To evaluate whether climatic stability influenced the genetic diversity, we investigated whether the genetic diversity of Cerrado sample locations was higher in the Quaternary stable areas compared to unstable areas. We classified all locations as belonging to stable or unstable areas based on the stability map. Comparisons were done for each locus using the non-parametric Kruskal-Wallis test in PAST, in which the stability was considered the explanatory variable and the genetic indices (percentage of nucleotide diversity or haplotype diversity) the response variables.
We performed Tajima´s D [66] and Fu´s Fs [67] neutrality tests and estimated the R2 statistic [68] to infer past population expansion considering the structure found by BAPS. All tests were performed in DnaSP using 10,000 coalescent simulations and a significance level of 0.05. We also inferred the dynamic of population size as a function of time performing a coalescent extended Bayesian skyline plot (EBSP) analysis implemented in BEAST 1.8.4 [69]. We used the TN93+I substitution model for ND2 and HKY+I model for the FIB-5. Analyses were done implementing the lognormal relaxed clock model, since the null hypothesis of equal evolutionary rate throughout the tree was rejected for both markers. We set a mutation rate of 9.0x10 -9 per base pair per year for ND2 [70] and 3.6x10 -9 per base pair per year for FIB-5 [71]. We used a total of 400 million runs and saved parameters every 40 thousand runs and looked at the posterior effective sample size (ESS) in TRACER 1.6 [72] using a threshold of 200. Finally, we plotted the median of population size in log scale through time and displayed the 95% high posterior density intervals using R.
We evaluated the genetic differentiation (F ST ) and the mitochondrial haplotype sharing between pairs of populations to evaluate whether the proposed biogeographic corridors served as dispersal routes for N. fasciata, considering: 1) MACA and PNCG-as an indicator of historical connection throughout the Andes region and/or Madeira River; 2) MACA and PATO, AGCV, or EEAE as an indicator of historical connections across central Amazonia; and 3) MACA and URUC or PNCM as an indicator of historical connection along the Atlantic coast. Genetic differentiation was estimated using F-statistics as previously described.

Species distribution models
A total of 158 occurrence records of N. fasciata were used to construct the species distribution models (Fig 2D and S1 Table), with most of these records being distributed in Brazilian territory and few of them in Paraguay, Suriname, and Bolivia. The climatic niche model for the present showed a high predictive power, based on the high average AUC value (0.906, S2 Table).
Based on paleoclimatic models, suitable areas during the LIG were more fragmented and restricted to the southern part of the distribution, compared with the projection for the present time (Fig 2). The LGM was characterized by the expansion of suitable areas, mainly in the central and northern Cerrado distribution compared to the LIG. The Holocene seems to be the less favorable climatic period for the species, characterized by a remarkable reduction of suitable areas compared to the LGM. Suitable area patches concordant with the central Amazonia belt of low precipitation were observed during the LIG and Holocene, while small suitable areas concordant with the Andean slopes were found only during the LIG.
The distribution of the stable areas over the last 120 kyr is roughly concordant with some Brazilian plateaus: Parecis, Paraná-Guimarães, Serra da Canastra, Chapadão do São Francisco and Central (Fig 3B). All of these plateaus have characterized the landscape since the Paleogene (43-20 million years ago-myr) [73,74]. In the last 20 myr, these plateaus have been dissected and denuded as a result of temperature and precipitation increases [74,75], forming depression compartments between the relict relief. The duality between depressions and plateaus generates discontinuous or fragmented compartments. In Central Brazil, the relict relief with the greatest territorial extension is identified as Core (Fig 3). Depressions and valleys are located between Core compartments where Interior and Edge classes can be identified (Fig 3B).
Currently, patches are observed in the north of the Brazilian Cerrado where the plateaus are less frequent in this region. Considering the distribution of the stable areas, the sampling sites AGCV, EEAE, PGSV, PNEM, PNCG, and PNSC were assigned to stable areas, while URUC, PATO, and PNCM corresponded to areas without climatic suitability (or unstable areas) in at least one of the projections (Fig 3A).
The ND2 haplotype network was represented by one broadly distributed and common haplotype connected by a few mutational steps with several less frequent haplotypes in a star-like pattern (Fig 4). This pattern can be interpreted as a signal of recent population expansion. Moreover, the ND2 network also showed a haplotype group composed exclusively of individuals sampled in PGSV and BRMI, and a second group comprised of individuals sampled in PNSC and PNEM. Both haplotype groups were separated from the common haplotype by several mutational steps. The FIB-5 haplotype network showed three of the most common haplotypes connected by a maximum of three mutational steps from each other and widely distributed across sample locations (Fig 4). The genetic structure obtained by BAPS analysis was similar for both markers and showed support for only one genetic group.
The results of Mantel tests for both ND2 (r = 0.0059, Z = 53.8917, p = 0.5218) and FIB-5 (r = -0.069, Z = 3.7704, p = 0.662) indicated that the geographic distances do not correlate with the genetic differentiation among sample locations. The analyses of the influence of geographic gradients in genetic diversity using the ND2 indicated a negative correlation between nucleotide diversity and latitude (rs = -0.728, p = 0.007), characterized by a south-to-north trend of declining genetic diversity (Table 2). For FIB-5 there is an east-to-west trend of declining nucleotide diversity (rs = 0.580, p = 0.048) ( Table 2). There is a significant difference between ND2 nucleotide diversity sample medians of stable and unstable areas, with greater values in stable areas compared to unstable ones (χ2 = 5.4, p = 0.020) ( Table 3). For the remaining comparisons, the differences between stable and unstable areas were not statistically significant (Table 3).
Based on the genetic structure found, neutrality tests were carried out considering all sample sites as a single genetic population. Deviations from neutrality were found in Tajima

Median joining networks and haplotype distribution maps obtained using NADH dehydrogenase subunit 2 gene (ND2) and beta-fibrinogen gene intron 5 (FIB-5) datasets.
In the networks, each haplotype is represented by a numbered and colored circle: colors correspond to the sample sites, and circle size and pie charts are proportional to the number of individuals (see Table 1 and Fig 1 to interpret). Mutational steps are represented by hatch marks and missing haplotypes are indicated by small black circles. In the haplotype maps, colors correspond to haplotypes in the sample site, and the size of circle and pie charts are proportional to the number of individuals (see Table 1 and Fig 1 to interpret). Non-shared haplotypes (unique) are shown in gray. https://doi.org/10.1371/journal.pone.0212876.g004

Marker
Genetic

Discussion
In the present study, we assessed the levels of genetic diversity and structure in Neothraupis fasciata based on samples from a wide geographic area in the Cerrado. Our data support a single population (supported by BAPS results) with a weak genetic structure (global F ST ), in which southernmost locations showed more genetic variability than those sampled in the northern Cerrado. Considering that areas of greater genetic diversity are commonly associated with the origin center of a species (e.g. [25,27,76]), it is plausible to assume that the origin center of the N. fasciata lineage is in the southern part of the Cerrado. The southern Cerrado encompasses some topographical relief core areas, such as Parecis, Paraná-Guimarães and Serra da Canastra plateaus, which are, at least in part, concordant with some predicted stable areas (see Fig 3). These historically stable areas may have played an important role in the speciation process and also served as refugia favoring higher levels of genetic diversity for this species. Additionally, according to the phylogeny of the Thraupidae family [77], the distribution range of the N. fasciata sister species (Lophospingus pusillus, Lophospingus griseocristatus, Gubernatrix cristata, and Diuca diuca) is restricted to southern South America [78]. The observed genetic structure in N. fasciata does not seem to be related to the isolationby-distance hypothesis. This is an unexpected result, given the low dispersal from the birth area reported for this species (around 300 meters; [37]). However, a previous study of the genetic structure of this tanager using microsatellite markers also found no effect of geographic distance under the genetic differentiation among these same sampled locations [39]. The isolation-by-distance pattern occurs when the migration-drift equilibrium is reached, but it can take a relatively long time [79]. Therefore, the demographic history of a species should be considered, since the isolation-by-distance pattern is most commonly found in older populations and may be absent in populations established in a more recent time [80]. Considering that our results suggest a recent population expansion in N. fasciata, it seems reasonable to assume that the mutation-drift equilibrium has not been reached yet, and this may be masking the effect of the low dispersion capacity of this species in the spatial distribution of its genetic diversity.
The Cerrado landscape is characterized by large blocks of plateaus (Fig 3) separated from each other by peripheral depressions, probably formed at least 2-3 million years ago [13]. Our ecological niche models associated with the topographical information showed that the major stable areas were associated with different plateau areas, akin to previous studies that found stable areas along central Brazilian plateaus [3,17,19,[81][82][83]. Although Ab'Saber [81] suggested the existence of a large and continuous stable area along the plateaus in the southern Cerrado during the Quaternary, this hypothesis was not supported by our findings (Fig 3) and by other Phylogeography of the White-banded Tanager studies that found small and isolated stable areas along the southern Cerrado [3,82,83], plausibly representing isolated species refugia. Our paleodistribution models also showed an increase in suitable area and connection from LIG to LGM, and a decrease and fragmentation from the LGM to the Holocene, followed by an expansion from the Holocene to the present distribution of N. fasciata. The shifts in the distribution of habitat suitability for N. fasciata presumably accompanied the shifts in the distribution of the Cerrado itself, since this species is strictly associated with cerrado sensu stricto and grassy cerrado [32][33][34][35]. In this sense, our results are consistent with the predicted expansion of South American open biomes during cold and dry periods, and retraction during warmer periods (e.g. [1,7,8,10,11,84]).
According to our results, N. fasciata underwent a recent population expansion event starting about 60 kyr (neutrality tests and Figs 4 and 5) likely favored by the increasing in area and connection of suitable habitats between LIG and LGM. Although the increase in connection of suitable areas was detected in the projection for 21 kyr, it is possible that this change started earlier. The period from 75 to 12 kyr (marine isotope stages 4-2) was marked by climatic changes reflecting a decrease of the average temperature, characterized by minor glacial cycles [85,86], which may have contributed to the increase in the connection of suitable habitats for N. fasciata. Although our results indicated an event of population expansion, studies encompassing open-vegetation vertebrates of South America showed different trends for the LGM, with events reported of decrease (mammals and birds), expansion (mammals), or even constant effective population size (mammals and lizards) (see [4]).
We found greater genetic diversity in those individuals from the southern sample sites: (I) PNCG, (II) PNSC + PNEM, and (III) PGSV + BRMI, and a low number sharing haplotypes with each other (Fig 4). These results suggest that these populations were partially isolated before the expansion event. Events of population expansion may result in a gradual reduction in genetic diversity toward newly colonized areas due to successive founder effects (leadingedge expansion model) and the increase in genetic diversity in populations that received Phylogeography of the White-banded Tanager migrants [26,27]. Thus, our results indicate that the process of population expansion promoted the reestablishment of gene flow in the southern part of the species distribution, reconnecting genetic groups historically isolated in different plateaus (Fig 3). Additionally, the gradual reduction in genetic diversity following a south-to-north trend indicates that the species may have survived in isolated refugia in the southern part of its distribution during less favorable periods with subsequent colonization events towards the northern Cerrado. This process was likely favored by the continuity of vegetation formed in the core Central plateau and a possible extension of the cerrado sensu stricto-like vegetation toward northern residual relief (see transitional and patch areas in Fig 3B) in the interplateau depressions promoted by the decrease in temperature between LIG and LGM. Plateaus have deep and well-drained soils that are favorable to the establishment of the cerrado sensu stricto vegetation, while interplateau depressions have younger soils and, consequently, a more heterogeneous landscape than the plateaus (see Silva 1995 for a review). As in the plateaus, intraplateau depressions are composed of planning surfaces covered predominantly by Oxisols, which favor the establishment of cerrado sensu stricto-like vegetation in response to climatic changes. In this regard, the interplateau depressions are possibly more prone to vegetation range shifts due to climatic changes than the plateaus, as denoted by our paleomodelling results and the literature that have associated the Cerrado stable areas with plateaus and residual relief [3,17,19,[81][82][83].
We observed higher genetic diversity and a signal of historical differentiation between populations located in different stable areas. Thus, the mitochondrial results reinforce the hypothesis that high levels of divergence and population structure are expected across refugia [31] and the genetic diversity is expected to be higher in stable areas compared to unstable areas [19,30,31,83]. In addition, we stress that N. fasciata sample locations categorized as unstable are in the northern part of the distribution, and the effects of range expansion may act as a confounding factor. Therefore, additional sampling mainly in unstable areas is necessary to disentangle the effects of population range expansion and historical stability.
Connections among the South American savannas are postulated to have occurred through biotic corridors: along the Andes region [20]; right across central Amazonia, following a belt of low precipitation [15]; and along the Atlantic coast [84] (see [1] for a review). Additionally, a fourth corridor crossing Amazonia from the southwest to the northeast following the Madeira River was postulated as an alternative dispersal route [18]. Our results showed that MACA and PATO hold the lowest level of genetic differentiation (S4 Table) and are the only ones to share haplotypes, suggesting that N. fasciata dispersal may have occurred mainly through central Amazonia corridor, besides we were not able to dismiss the use of the other hypothesized corridors. Additionally, the niche models indicated a region of habitat suitability concordant with the low precipitation belt during the LIG. In accordance, Bueno et al. [17] and Buzatti et al. [19] reported the existence of suitable climatic conditions to establish different Cerrado tree species throughout central Amazonia, reinforcing the existence of the central Amazonian corridor. However, our results are not in accordance with the literature that suggests the corridors along the Atlantic coast [1] or the Madeira River [18] as the main dispersal routes used by N. fasciata. The pattern observed in our results and the growing literature about the importance of the central Amazonian corridor (e.g. [17,19]) highlight the need for more studies to understand the importance of the hypothesized Quaternary savanna corridors to the South American avifauna dispersal.

Conclusion
Genetic groups of N. fasciata located in the southern part of the species distribution (likely its center of origin) might have persisted in partially isolated populations restricted to the plateaus during the cycles of savanna contractions. Subsequent events of savanna expansion promoted the reestablishment of gene flow among these isolated populations and the increase in the distribution range through the (re)colonization of the northern portion of the species distribution. Thus, the intrinsic relationship between the relief heterogeneity (plateaus and depressions) and the climatic fluctuations, mainly in the Pleistocene, promoted population reconnection and demographic expansion of N. fasciata, reinforcing previous studies (e.g. [82]) that suggest that this relation of cause and effect (the interaction between climate and relief defining vegetation dynamic) may have been decisive in the diversification of the Cerrado biodiversity.
Supporting information S1