Evolutionarily significant units of the critically endangered leaf frog Pithecopus ayeaye (Anura, Phyllomedusidae) are not effectively preserved by the Brazilian protected areas network

Abstract Protected areas (PAs) are essential for biodiversity conservation, but their coverage is considered inefficient for the preservation of all species. Many species are subdivided into evolutionarily significant units (ESUs) and the effectiveness of PAs in protecting them needs to be investigated. We evaluated the usefulness of the Brazilian PAs network in protecting ESUs of the critically endangered Pithecopus ayeaye through ongoing climate change. This species occurs in a threatened mountaintop ecosystem known as campos rupestres. We used multilocus DNA sequences to delimit geographic clusters, which were further validated as ESUs with a coalescent approach. Ecological niche modeling was used to estimate spatial changes in ESUs’ potential distributions, and a gap analysis was carried out to evaluate the effectiveness of the Brazilian PAs network to protect P. ayeaye in the face of climate changes. We tested the niche overlap between ESUs to gain insights for potential management alternatives for the species. Pithecopus ayeaye contains at least three ESUs isolated in distinct mountain regions, and one of them is not protected by any PA. There are no climatic niche differences between the units, and only 4% of the suitable potential area of the species is protected in present and future projections. The current PAs are not effective in preserving the intraspecific diversity of P. ayeaye in its present and future range distributions. The genetic structure of P. ayeaye could represent a typical pattern in campos rupestres endemics, which should be considered for evaluating its conservation status.

. From a population viewpoint, a species is not always a single unit for conservation because intraspecific subunits may evolve independently in response to climate disruptions (Bálint et al., 2011;Forester, DeChaine, & Bunn, 2013). Due to these reasons, current assessments about the future effectiveness of PAs could be over-optimistic (Pauls et al., 2013). Therefore, an effective plan for future species viability should incorporate the intraspecific levels of biodiversity. In this context, evolutionarily significant units (ESUs) represent ideal targets for conservation because they contain the raw material for future evolutionary adaptations (Bowen & Roman, 2005;Fraser & Bernatchez, 2001;Pauls et al., 2013). Thus, population viability evaluations should begin with the challenging task of delineating these intraspecific units (Fraser & Bernatchez, 2001).
The conservation status of amphibians is alarming in this context because they are among the most threatened vertebrates (Pimm et al., 2014) and many species are suffering population declines associated with a variety of threats, including climate change (Blaustein et al., 2011;Stuart et al., 2004). Despite these serious conservation concerns, 42% of the amphibian richness is misrepresented or completely outside PAs (Nori et al., 2015). Amphibians generally exhibit high levels of genetic structuring usually associated with the geographic isolation of demes (Allentoft & O'Brien, 2010;Rodríguez et al., 2015), which makes it difficult to assess whether the differentiated lineages are either distinct species or populations (e.g., Carnaval & Bates, 2007;Gehara, Canedo, Haddad, & Vences, 2013). In these cases, if conservation policies do not consider lineage differentiation, a significant loss of cryptic diversity is expected under climate change (Bálint et al., 2011;Pauls et al., 2013). For this reason, understanding intraspecific ESUs, their evolutionary history, spatial distribution, and connectivity are important steps toward minimizing genetic diversity loss among isolated amphibian populations (Bálint et al., 2011;Beebee & Griffiths, 2005;Pauls et al., 2013;Weeks, Stoklosa, & Hoffmann, 2016).
The reticulated leaf frog Pithecopus ayeaye B. Lutz, 1966 (Anura, Phyllomedusidae) (Figure 1) is a high-altitude endemics leaf frog of southeastern Brazil, which reproduces throughout the rainy season (from October to January) (Oliveira, 2017). The species is rare and shows low individual density in most of its range (Araujo, Condez, & Haddad, 2007;Baêta, Caramaschi, Cruz, & Pombal, 2009). It produces clutches with few large eggs and the larvae generally inhabit rocky stream pools, with crystal clear and slow-flowing water (Oliveira, 2017;Pezzuti, Leite, & Nomura, 2009). The males are territorials, defending pulleys along streams with riparian vegetation for oviposition where funnel nests are constructed above the water using the leaves (Nali, Borges, & Prado, 2015;Oliveira, 2017). These leaves are not randomly chosen, but those from Melastomataceae bushes on the margins of the rivulets are preferred (Oliveira, 2017), probably due to foliar area that minimizes the risk of clutch desiccation and the presence of spiny  (Dias, Maragno, Prado, & Cechin, 2014;Oliveira, 2017). When tadpoles hatch from eggs, they fall into the water and complete their development (Dias et al., 2014).
The species is susceptible to stochastic climatic events, like dry periods in the rainy season that can cause severe stream drought and dehydration of egg clutches.
Pithecopus ayeaye is classified as critically endangered (CR) by the International Union for Conservation of Nature (IUCN; Caramaschi, Cruz, Lima, & Brandão, 2010) based on B1ab(iii)+2ab(iii) criteria, as it was only known from two disjunct localities threatened by severe habitat loss due to mining and human-induced fires (Caramaschi et al., 2010). The threat sources are more numerous because the mountain grasslands where P. ayeaye occurs are also menaced by forestry, cattle farming, nonsustainable "ecotourism," and poorly planned urbanization (Silveira et al., 2016). These threats impact directly or indirectly on the highland streams where P. ayeaye breeds. Hence, P. ayeaye was included as a priority in a Brazilian national action plan for Espinhaço Range herpetofauna (ICMBio, 2011). After the discovery of new localities (Araujo et al., 2007;Baêta et al., 2009), the species no longer meets the geographic requirements to be categorized as CR. For this reason, P. ayeaye was recently removed from all threat categories in the Brazilian List of Endangered Species (ICMBio, 2014), and excluded from the national agenda of conservation priorities. However, the impacts of the fragmented distribution of the species on its genetic diversity and adaptive potential were not considered in the current categorization.
Pithecopus ayeaye occurs in a threatened ecosystem (Fernandes, Barbosa, Negreiros, & Paglia, 2014) and exhibits a naturally fragmented distribution in mountaintops of southeastern Brazil. Although IUCN informs that the species occurs at elevations over 647 m a.s.l. (Caramaschi et al., 2010), probably due to a misunderstanding of Araujo et al. (2007)'s description of Furnas do Bom Jesus State Park limits, literature and collection data suggest that the distribution is mainly restricted at elevations higher than 900 m a.s.l. (Table S1). It occurs mainly in the southern limits of the campos rupestres grasslands (see Silveira et al., 2016), but it is also found in grassland patches of the Poços de Caldas Plateau (Araujo et al., 2007;Baêta et al., 2009).
Nevertheless, it has been largely neglected in Brazilian research and conservation policies (Jacobi et al., 2007;Silveira et al., 2016), even though there is clear evidence of an ominous future for this unique mountaintop ecosystem. For example, a climatic model has predicted a loss of 95% of suitable areas, in an optimistic scenario, by the year 2080 (Fernandes et al., 2014).
This intraspecific diversification has been frequently used to identify conservation units within species, which have been applied in the definition of Evolutionarily Significant Units (ESUs). Moritz (1994) used phylogeographical analysis based on sequence markers, defining defined ESUs as reciprocally monophyletic, mitochondrial DNA (mtDNA) groups with significant divergence in nuclear alleles frequencies among them. The advantages of this criterion are its objectivity and generality, applicable to phylogeographical data (Fraser & Bernatchez, 2001), but it might not be able to distinguish ESUs with incomplete lineage sorting (ILS) or mtDNA gene flow. To accommodate these genealogical processes, Fraser and Bernatchez (2001) proposed that ESUs are intraspecific lineages with highly restricted gene flow among them, allowing ESUs delimitation without reciprocal monophyly.
In this study, we aimed to delimit ESUs within P. ayeaye using a broad geographic sampling and multilocus sequence data in a statistical phylogeography framework coupled with GIS information (e.g., Forester et al., 2013). We analyzed multiple island models to validate ESUs, taking ILS and gene flow into account, and used niche divergence tests to refine predictions about the effectiveness of PAs in preserving diversity at the present and in the future (Beerli & Palczewski, 2010;Carstens, Brennan, et al., 2013;Warren, Glor, & Turelli, 2008).
In order to evaluate the future survival of distinct ESUs in the current PAs, we used ecological niche modeling (ENM) to project ESUs' ranges based on the ongoing and future scenarios of climate change (Bálint et al., 2011), applying rigorous model evaluation to deal with the scarcity of occurrence data due to the rarity of the species (Shcheglovitova & Anderson, 2013). From the available evidence, we expect a spatially associated intraspecific diversification for P. ayeaye and the loss of suitable areas for the species is expected in a future scenario of climate change. Here we assume that the geographic distribution of the species is partially explained by climate, as the areas where P. ayeaye occurs are islands of subtropical climate in a tropical matrix (Alvares et al., 2013), which could favor a set of physiological and/or behavioral adaptations.

| Data collection
We obtained tissue samples mainly from museums and collections, many of them with voucher specimens (see Appendix S1). These were complemented with samples collected by us from a few other locations to bring our total to 88 individuals from 13 sample points in nine counties ( Figure 2, Table S1), which is a significant sampling in relation to known species geographic distribution (see Baêta et al., 2009). In the case of tadpole sampling, we excluded individuals in similar developmental stages collected together in the same place or stored in the same museum lot, in order to avoid first-degree relatives and, consequently, to maximize sampling randomness. We extracted genomic DNA from liver or muscle samples using the phenol-chloroform protocol (Sambrook & Russel, 2001). For all specimens, we amplified and RPL3intR, . Polymerase chain reactions (PCRs) were performed in a 15 μl reaction volume containing: 20 ng of genomic DNA, 1× buffer, 2.5 mmol/l MgCl 2 , 1.25 μmol/l each primer, 3 mmol/l dNTPs, 0.72 μg bovine serum albumin, and 0.625 U Platinum™ Taq DNA polymerase (Thermo Fisher Scientific). Amplifications were performed as one initial denaturation for 95°C during 5 min followed by 35 cycles [denaturation at 95°C for 30 s, variable melting temperatures and times between fragments (54°C by 40 s to Cyt-b; 60°C by 50 s to POMC; and 62°C by 40 s to RPL3), extension at 72°C for 1 min/1,000 bp] and a final extension at 72°C for 7 min. We reduced temperature melting, increased time melting, and/or increased up to 3.5 mmol/l MgCl 2 for samples difficult to amplify. All single PCR products were purified using a Sambrook and Russel's (2001) polyethylene glycol 20% protocol, with some modifications (Santos Júnior, Santos, & Silveira, 2015).
Purified amplicons were fluorescence-marked through BigDye™ TerminaTor v3.1 cycle sequencing kit (Thermo Fisher Scientific) following the manufacturer's instructions. The primers used in this reaction were the same to those used in PCRs, except for RPL3-intF , which was replaced with RPL3-P3 (5′ WCTGGCCTGCTCTGGTTAT 3′) designed by us in Primer3Plus (Untergasser et al., 2007). The marked amplicons were analyzed in an automatized ABI 3130xl DNA sequencer (Thermo Fisher Scientific) in both directions.

| Sequence edition and data characterization
Sequence fluorograms were interpreted, assembled, pre-aligned, and edited in SeqScape™ 2.6 (Thermo Fisher Scientific). Edited fragments were aligned with the cluSTalW module of the mega7 software (Kumar, Stecher, & Tamura, 2016;Larkin et al., 2007), with gap opening penalized 10 times more than extension for intron alignment. The gametic phases of nuclear markers were reconstructed through the algorithm implemented in phaSe 2.1.1 software (Stephens, Smith, & Donelly, 2001), which were interconverted to fasta format using SeqPHASE web tool (Flot 2010). In cases where more than one haplotype pair was reconstructed, we selected the more probable pair for subsequent analyses, except those based in genealogies reconstruction. In these cases, haplotype phases not fully resolved (PP < .9) were checked by eye and the uncertain nucleotide positions were kept. Because there were some heterozygous individuals for insertions and deletions in RPL3, we estimated their haplotypic phases via the codification of superimposed traces in inDelligenT 1.2 web-based software (Dmitriev & Rakitov, 2008).
After these steps, we selected the best fit model of molecular evolution among 88 substitution schemes using the BIC criterion in jmoD-elTeST 2.1.10 software (Darriba, Taboada, Doallo, & Posada, 2012) for each DNA fragment.
We did a maximum likelihood (ML) test of the strict molecular clock for each DNA fragment under the same substitution model selected in the previous step as implemented in mega7 (Kumar et al., 2016). Because an input tree is required for this test, we generated gene trees of unique haplotypes in RaxML 8.2.4 software (Stamatakis, 2014), choosing the "best scoring ML-tree with rapid bootstrap" option under the GTRCAT model. The null hypothesis of equal evolutionary rate along the tree was not rejected for both Cyt-b (p = .85) and POMC (p = .47), but it was rejected for RPL3 (p < .001). We applied this test to build simpler and less parametric models in subsequent analyses, aiming to achieve a balance between bias and variance in results (Kelchner & Thomas, 2006).

| ESUs discovery, diversity, and relationships
We applied a two-step ESU delimitation consisting of sequential stages of discovery and validation, which is analogous to the rationale proposed by Carstens, Pelletier, Reid, and Satler (2013) for species delimitation, but accommodating gene flow explicitly. For the ESUs discovery, we first tested spatial population structure using the genelanD package in r (Guillot, Estoup, Mortier, & Cosson, 2005;Guillot, Mortier, & Estoup, 2005; R Core Team 2016). We chose a spatially explicit model because of our expectation of structure is related to isolation in sky islands. The mtDNA data were entered as haplotypes, while the nDNA haplotypes were encoded as alleles for the analysis. Because our goal was to identify those lineages that minimized global gene flow (Fraser & Bernatchez, 2001), we opted for the uncorrelated allele frequencies model, which is generally unable to detect subtle structuring (Guillot, 2008). We made ten parallel runs with 1 × 10 6 iterations each. The analysis iterations each, with a thinning of 1 × 10 3 , and K varying between 1 and 10 biogeographical units (BUs), which we consider as being putative ESUs. Additionally, we assessed the isolation by distance (IBD) hypothesis to verify whether genetic diversity can be explained by geographic distance between localities. For this, we performed a Mantel test between log 10 -transformed geographic and genetic distances matrices of our 13 sample points in iBDwS (Jensen, Bohonak, & Kelley, 2005), with 10,000 randomizations. The multigenic distance matrix was calculated in pofaD (Joly & Bruneau, 2006) using the genpofad algorithm (Joly, Bryant, & Lockhart, 2015).
To estimate the relationship among discovered BUs, we constructed a lineage tree (=species tree) in the STarBeaST2 (Bouckaert et al., 2014;Ogilvie, Bouckaert & Drummond, in press) using a Yule prior with a constant population model. To estimate divergence times, we used an uncorrelated log-linear clock model for RPL3, and strict clock models for POMC and Cyt-b. In the latter case, we used a standard mtDNA substitution rate (mean of 0.01 substitutions per lineage per million years; Johns & Avise, 1998), due to the lack of fossils for calibration. This is similar to the ND2 rate (Crawford 2003), which is widely used in amphibian dating (e.g., Carnaval & Bates, 2007). This analysis was made with two replicates, with a pre-burn-in of 2.5 × 10 7 followed by 7.5 × 10 7 iterations each. The analysis was repeated twice, and we checked for convergence, stationarity, and minimum adequate effective sampling size (ESS > 200) of analyses in Tracer v1.6 (Rambaut, Suchard, Xie, & Drummond, 2014).
Furthermore, we built statistical parsimony networks in popART 1.7 software (Leigh & Bryant, 2015;Templeton & Sing, 1993) to visualize the relationships between the haplotypes of each gene fragment.

| ESUs validation
We estimated the historical gene flow (effective number of migrants, M = 4N e m) and the genetic diversity (ϴ = 4N e μ) of the populations under a coalescent framework using migraTen (Beerli, 2006). We first implemented a Bayesian full model accounting for the BUs (Fig. S1, model 1) using an empirical transition-transversion ratio (R Cyt-b = 4.661; R RPL3 = 1.154; R POMC = 1.449) estimated under the K80 model in mega7 (Kimura, 1980;Kumar et al., 2016), and a scheme of relative mutation rates among loci. Individuals with missing data in haplotypes were excluded, as this may result in spurious results in migraTen (Carstens, Brennan et al., 2013). Initial values for parameters were derived from F ST estimates. We conducted the analysis with two parallel runs, using a static heat strategy, setting one long and 12 short chains, with the cold chain equal to one, the hottest chain equal to 5 × 10 5 and values of the remaining chains growing at a cumulative exponential scale of x 1.4 , starting from 1.5.
We made a previous burn-in of 1.25 x 10 6 generations, and 2.5 × 10 6 states were visited in each run, with a thinning of 100. Exploratory analyses with the full model allowed us to determine the sampling window for ϴ and M, and we assumed a normal distribution for these parameters.
Based on the genelanD and STarBeaST2 results, we generated a set of seven reduced models to validate the putative ESUs assignment ( Fig. S1, models 2-8). Fraser and Bernatchez (2001) advocated the use of criteria that "provide evidence of lineage sorting through highly reduced gene flow," but they did not propose a cutoff level of gene flow for this criterion. Therefore, we tested scenarios of BU independence against scenarios of split between sister BUs, and split of all BUs (Fig. S1). In summary, we tested models ranging between all BUs forming a single ESU, and each BU as an independent ESU, with distinct routes of gene flow among them. If the gene flow is so high that two or more BUs form a single genealogically cohesive cluster, the BUs are collapsed into a single ESUs. In order to reduce the potential set of hypotheses, they were constructed under the following assumptions: (1) when present, gene flow is always bidirectional, and (2) sister BUs have gene flow between them, except in isolation models (see Fig. S1 for a graphical representation). The second assumption was based on the expectation that gene flow tends to decrease with longer splitting times between lineages (Pinho & Hey, 2010). Besides that, we take into account the nearest shared node from the lineage tree to collapse the BUs. We calculated the marginal likelihood of each model and compared them under a Bayes factor test using Bezier's approximation score (Beerli & Palczewski, 2010). Under this approach, we can select the model with the highest probability of fitting our data. We used model averaging for the final estimates of ϴ and M, and evaluated the performance of MCMC sampling through ESS and acceptance ratio values.

| Parameters estimation through Approximate Bayesian Computation (ABC)
Because STarBeaST2 analyses do not take into account potential gene flow, and consequently, it can underestimate divergence time (τ) among ESU (Pinho & Hey, 2010), we also implemented an Approximate Bayesian Computation (ABC) approach to co-estimate demographic parameters based on the lineage tree topology and the gene flow routes obtained from the best island model of ESUs.
The older the lineage divergence, greater lineage sorting is expected (Maddison, 1997). For this reason, the aim of the tree dating analysis was to verify the congruence between the selected models and the time of ESU's diversification. We simulated 1 × 10 6 coalescent genealogies for three loci with mS (Hudson, 2002), and processed them to obtain the following summary statistics in the mSSS.pl script (Takebayashi, 2011): average nucleotide pairwise distances per locus (π), number of segregating sites (S), Tajima's D (Tajima, 1989), π within populations, and π between populations. We used uniform prior distributions for all parameters including ϴ per locus (lower bound: 0.01, upper bound: 10), divergence times in 4N e units (0.0001, 0.5), and migration rates in 4N e m units (0, 50). The same summary statistics were calculated for the empirical data globally, and for each ESU and locus in DnaSP (Table 1, Table S2). DNA divergence between ESUs (π between populations) was calculated as the average number of nucleotide substitutions per site between populations (Nei, 1987) for each locus in DnaSP (Librado & Rozas, 2009) (Table S2)

| Ecological Niche Modeling (ENMs)
Using the bioclimatic data at a 2.5-min resolution of latitude and longitude, we modeled species niche and projected distribution of P. ayeaye across the mountaintops of a selected region in southeastern Brazil. Because we had many more collection records than sample points (Table S1), we used the putative population limits estimated by genelanD to assign known records to our delimited ESUs (Fig. S2). We obtained the climate layers from Hijmans, Cameron, Parra, Jones, and Jarvis (2005) to characterize the environmental space for ENMs using both present  and future cli-  (Table S2). This method is based on the correlation matrix among variables to minimize collinearity problems, consequently avoiding biased predictions of the ENMs.
A key assumption in our modeling is that ensemble forecasts based on multiple models generate more accurate or at least more conservative projections of species distribution than single models do (Araújo & New, 2007). Therefore, four modeling methods were used to build the ENMs, including Bioclim, Gower Distances, Maximum Entropy (MaxEnt), and Support Vector Machine (SVM), both presence-only or presence background methods, and adequate to our limited data by its simplicity in the configurations (for a review, see Peterson et al., 2011). All ENMs were implemented in the DiSmo package in r (available at https://CRAN.R-project.org/ package=dismo). In summary, models were first generated for present climate and then projected onto future conditions to predict the species geographic range for these two distinct time periods. We assessed model performance for each decision threshold using the "leave-one-out test" because of the small number of occurrence records for P. ayeaye and ESUs (Table S1). Hence, multiple predictions were made for P. ayeaye and ESUs, selecting one occurrence record for removal in each case. This approach is described as a variation to the k-fold partitioning method on which a Jackknife sampling is imposed (Bean, Stafford, & Brashares, 2012;Pearson, Raxworthy, Nakamura, & Peterson, 2007;Shcheglovitova & Anderson, 2013).
For each prediction, we applied the lowest presence decision threshold (LPT) to test the ability to predict the deleted occurrences. If the ENM successfully predicts both a small area and the deleted occurrence record, it is better than a random model (p < .05). However, if the model predicts a large area and fails to predict the deleted occurrence record, it is not considered a good model (p > .05). Therefore, the p-value is calculated from the success and failure ratio of the prediction (Pearson et al., 2007).
Our modeling procedure resulted in combined suitability maps of the best models for P. ayeaye and ESUs, and for each climate condition (4 ENMs* occurrence records, for present climate and 4 ENMs * 3 AOGCMs for future climate condition * occurrence records, for future conditions). Finally, we obtained a consensus map for each combination of ENMs and AOGCM. To assess individual model variability sources, we separated and mapped uncertainties in forecast ensembles (Diniz-Filho et al., 2009). For this purpose, we performed a twoway ANOVA for each grid cell using suitability as response variable and the methodological components (AOGCMs and ENMs) as explanatory variables.

| Effectiveness of PAs to ESUs protection
We were interested in verifying the PAs effectiveness to protect the entire P. ayeaye's suitable area and each ESU separately. We over- T A B L E 1 Global and population summary statistics for sampled loci. Parameters shown are the total number of haplotypes (N), number of unique haplotypes (h), haplotype diversity (Hd), and nucleotide diversity per site (π). Tajima's D values were explained by chance

| ESUs characterization and relationships
We found correlation between genetic structure in P. ayeaye and mountain ranges where it occurs. Geneland runs showed convergence in posterior probabilities (PP) of models, after a burn-in of 1 × 10 3 replicates. The analyses returned three BUs with PP = .47 (Fig. S3).
The geographic distribution of these lineages corresponded to the 2) had PP = .33 (Table 2). This non-negligible probability for an alternative model may be due to incongruence among loci. While Cyt-b and RPL3 show 76 and 96% of probability associated with model 3, respectively (Fig. S1), the best fit model for POMC was number 2, with 43% of probability (Fig. S1). In any case, our results suggest higher gene flow between the sisters ESUs of Canastra and Poços, but limited gene flow between Quadrilátero and at least one of the other ESUs.
Furthermore, models of island isolation and panmixia had no support in either single-or multilocus analyses (Table 2). From the genetic point of view, the three BUs are ESUs from a historical metapopulation with a stepping-stone island pattern among them.
Regarding ABC estimates, the evaluation of model fit done with a PCA analysis of summary statistics suggested that the simulated demographic model fits well the empirical data because the cloud of posterior data points matches the observed data point (Fig. S4). As expected, the divergence times estimated with the ABC approach (which takes gene flow among ESUs into account) were older than the STarBeaST2's estimates, especially the older splitting event, which was about 3.5 times larger on average (Table 3). Furthermore, ϴ estimates had values within the 95% highest posterior density (HPD) confidence intervals of migraTen parameters, except for Canastra's ϴ (Table 3). On the other hand, all migration rates were higher than those estimated in the full Bayesian approach (Table 3) The Canastra ESU exhibited much more genetic diversity than the others (Tables 1 and 3). This does not seem to be a sampling artifact effort and geographic coverage (Table S1, Figure 3). The RPL3 marker presented excess of haplotype pairwise differences due to divergence between geographically divided ESUs (Figure 4c), further strengthening our BUs as ESUs (Tables 1 and 3).

| ENMs and effectiveness of PAs to ESUs protection
We used the four ENM methods to build a consensus map of the projected present and future potential distribution for entire P. ayeaye, Canastra and Quadrilátero ESUs ( Figure 5). Poços ESU was not included due to the scarcity of occurrence data to construct reliable models (Table S1). The projections presented for P. ayeaye and BUs were trained using between seven to 26 localities and show high success rates in Jackknife tests (Table 4) The two-way ANOVA applied shows that the median of the variation projected in future climate conditions for P. ayeaye and Quadrilatero ESU are due to differences in ENMs (up to 85.6%; Table 5), with the lowest differences in the southern limit of Minas Gerais State (Fig. S4). Furthermore, the maps showing this sum of squares ( Fig. S5) indicate the largest differences among methods for P. ayeaye, Canastra and Quadrilatero ESUs. When considering the ENMs, both P. ayeaye (all populations) and individual ESUs can be potentially found in other PAs, in addition to those already known to protect the species.
However, <15% of all suitability areas (even in more relaxed cutoff) are in protected areas for present projections of the species (Figure 6).
For the Canastra ESU, more than 35% of highly suitable areas are within PAs, but this is drastically reduced in future conditions (~4%) ( Figure 6). Even those models predict wide suitable climatic areas in PAs for the current climate, the known occurrence present only in rock outcrops should prevent its expansion in future climate conditions.

| DISCUSSION
Our results show that the P. ayeaye displays a genetic structure associated with the spatial fragmentation of a sky island environment.
Moreover, despite the reported wide range of P. ayeaye (Baêta et al., 2009) (Araujo et al., 2007;Baêta et al., 2009; FSFL, personal observation). Together, those five PAs encompass <4% of the entire, potential and suitable geographic distribution of P. ayeaye. In the future scenario, with the predicted expansion of the species, the coverage of protected areas for P. ayeaye tends to worsen or to improve unsatisfactorily, depending on suitability cutoff ( Figure 6). To aggravate the situation, P. ayeaye exhibits a set of life history traits associated with vulnerability in amphibians. Its rarity, K-reproductive strategy, and high habitat specialization are characteristics associated with low genetic diversity, increased population structure, and/or high threat levels (Cooper et al. 2008;Romiguier et al., 2014;Toledo, Becker, Haddad, & Zamudio, 2014;Rodríguez et al., 2015). Indeed, the smallest genetic diversity was observed in Poços and Quadrilátero ESUs, which may be associated with increased risks of local extinction. As P. ayeaye displays a naturally fragmented distribution of populations with evolutionary independence and is endemic to Brazil, we suggest that the conservation status of the species should be revised in the national list of endangered species because its long-term survival depends on local policies. For accurate recategorization, more research on demographic size and population trends needs to be conducted.
Our protocol to identify ESUs allowed us to validate P. ayeaye's ESUs even with the lack of reciprocal mtDNA monophyly (Moritz, 1994 (Carstens, Pelletier, et al., 2013). Geographically isolated populations may experience cyclical events of gene flow, especially in those species associated with interglacial refugia (Bonatelli et al., 2014). Gene flow breaks differentiation between lineages and generates mtDNA para-or polyphyly, a condition that does not meet the requirements of Moritz's (1994)  with the lack of well-defined haplogroups in the Cyt-b network. A varying level of haplotype or haplogroup sharing, as we observed between the ESUs, is expected as a result of historical gene flow before the more recent divergence event (Pinho & Hey, 2010). Therefore, the distribution pattern of P. ayeaye could be associated with interglacial refugia, a pattern that may be common in endemic species of campos rupestres (Bonatelli et al., 2014). This raises the hypothesis of potential gene exchange during historical phases of expansion and recontact between ESUs. Nevertheless, more detailed phylogeographical studies should be carried out to test this hypothesis, because knowledge about past distribution dynamics can provide insights for improved conservation decisions (Forester et al., 2013). Our results also suggest that the ESUs' delimitation method suggested by Moritz (1994) may be too restrictive, and more relaxed criteria (but not less methodologically rigorous) could be applied to phylogeographical data in order to accommodate cases of recent diversification and/or repeated recontact among formerly isolated lineages.
Despite this, the campos rupestres is expected to contract in warmer climatic conditions (Bonatelli et al., 2014;Fernandes et al., 2014), and our results show a surprising potential future expansion for the species' distribution. Thus, we can conclude that the major threats to the species are more related to habitat loss due to environmental degradation than to climate changes. The streams where the species breeds are extremely susceptible to erosion due to steep slopes and superficial soil. Currently, cattle grazing and degradation of marginal vegetation are the main threats of these fragile environments (e.g., Brandão & Álvares, 2009). In a more distal analysis, if climate changes lead to landscape structure alterations with negative impacts on stream rivulets, P. ayeaye could even be more threatened, even though our ENMs show an optimistic scenario for the species.
The niche overlap analysis indicates no climate selection between ESUs, which allows us to suggest conservation strategies that may be implemented in the present with potential positive impacts in the future. In situ conservation is the most feasible alternative for amphibian protection in Brazil (Haddad, 2008 (Corander, Waldmann, & Sillampaa, 2003). For example, concerning the Pedregulho and Alpinópolis sample points from the southern river-side of Rio Grande valley, there are only three samples and this large river may be a possible barrier to gene flow. Thus, we cannot rule out that increased sampling at this region might reveal another ESU. Even with these limitations, this study is important for providing valuable data about a species whose conservation status shows inconsistencies between global and national red lists. Finally, we used three independent sequence markers to make inferences about intraspecific diversification and relationships between P. ayeaye lineages. However, the use of a larger number of markers, like microsatellites and single nucleotide polymorphisms, should increase statistical power for clustering analysis and gene flow estimates (Allendorf, Honehlohe, & Luikart, 2010).
Our study sheds new light on conservation practices for amphibians in Brazil, as we found significant divergence among P. ayeaye populations that define three ESUs associated with distinct mountain regions, including one ESU found exclusively in an area without any kind of protection. This population structure reflects also the patterns found in many other taxa of the campos rupestres such as endangered lineages of the cactus Pilosocereus aurisetus (Bonatelli et al., 2014), the near threatened shrub Lychnophora ericoides (Collevatti et al., 2009), and the anurans Pithecopus megacephalus and Bokermannohyla saxicola (RFM, FRS and PCAG, unpublished data). Both have an associated structure with mountain ranges or geological subdivisions within them, as in P. ayeaye. Therefore, common preservation strategies can be applied for many Brazilian mountaintop endemic species to ensure their future viability. To accomplish this goal, intraspecific studies of the Brazilian sky islands' biota are greatly needed in order to guide decision-makers in generating policies that consider evolutionary and ecological processes and specificities of this ecosystem.

ACKNOWLEDGMENTS
We would like to give thanks to all funding agencies for financial support. Thanks to our laboratory colleagues for help in fieldwork,

CONFLICT OF INTEREST
None declared.

AUTHOR CONTRIBUTIONS
RFM, AC, PL and FRS designed the study; RFM, HT, RAB and FSFL collected samples; RFM, HT and FRS produced the molecular data; RFM, AC, PL and UO analyzed the data; RFM led and all authors contributed to the writing of the manuscript.

DATA ACCESSIBILITY
DNA data are archived in GenBank (accession numbers MF158349-158608). Information on the samples and occurrence points used in this study is available in Supporting Information.

SUPPORTING INFORMATION
Additional Supporting Information may be found online in the supporting information tab for this article.