Effects of climate changes on distribution of Eremanthus erythropappus and E. incanus (Asteraceae) in Brazil

Phylogeographic patterns of endemic species are critical keys to understand its adaptation to future climate change. Herein, based on chloroplast DNA, we analyzed the genetic diversity of two endemic and endangered tree species from the Brazilian savanna and Atlantic forest (Eremanthus erythropappus and Eremanthus incanus). We also applied the climate-based ecological niche modeling (ENM) to evaluate the impact of the Quaternary climate (last glacial maximum ~ 21 kyr BP (thousand years before present) and Mid-Holocene ~ 6 kyr BP) on the current haplotype distribution. Moreover, we modeled the potential effect of future climate change on the species distribution in 2070 for the most optimistic and pessimistic scenarios. One primer/enzyme combination (SFM/HinfI) revealed polymorphism with very low haplotype diversity, showing only three different haplotypes. The haplotype 1 has very low frequency and it was classified as the oldest, diverging from six mutations from the haplotypes 2 and 3. The E. erythropappus populations are structured and differ genetically according to the areas of occurrence. In general, the populations located in the north region are genetically different from those located in the center-south. No genetic structuring was observed for E. incanus. The ENM revealed a large distribution during the past and a severe decrease in geographic distribution of E. erythropappus and E. incanus from the LGM until present and predicts a drastic decline in suitable areas in the future. This reduction may homogenize the genetic diversity and compromise a relevant role of these species on infiltration of groundwater.


Introduction
Future climate change will create several impacts on biodiversity (Ravenscroft et al. 2015). A recent meta-analysis demonstrated that the negative impacts of habitat loss and fragmentation have been disproportionately severe in areas with high temperatures during the warmest month and declining rainfall (Zuber and Villamil 2016). Several studies addressing the impacts of climate change has been developed mainly focusing the loss of species diversity. Nevertheless, a few studies have evaluated the influence of climate on the genetic characteristics of populations.
Past bioclimatic scenarios models have shown that current species distribution follows a pattern of gene sharing since the Late Glacial Period (Collevatti et al. 2013). The cyclical climatic changes during the Quaternary period had a complex impact on the population dynamics Project funding: This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES)-Finance Code 001.
The online version is available at http://www.springerlink.com.
Electronic supplementary material The online version of this article (https://doi.org/10.1007/s11676-019-00968-z) contains supplementary material, which is available to authorized users.
& Lucas Fernandes Rocha lucasrochamg@gmail.com of tropical, Neosavanna species, consequently influencing their current species distribution and genetic diversity (Collevatti et al. 2015). While some species have shown resilience facing new climatic conditions (Kremer et al. 2012), future climate change will influence the species' genetic diversity and adaptation to new environmental conditions (Lima et al. 2017). The Brazilian savanna and the Atlantic forest are currently considered as global biodiversity hotspots in reason of the high level of susceptibility and the presence of endemic species (Myers et al. 2000). Additionally, Brazil's savanna is also classified as one of the most threatened savanna ecosystems on Earth (Gries et al. 2012) and the Atlantic forest biome is known as the hottest biodiversity hotspot, as a result of its high endemism of plants and vertebrates (Myers et al. 2000). Moreover, Brazilian high altitudinal grasslands have enormous numbers of endemism of plants, with restricted to isolated patches (Iganci et al. 2011). Besides, relevant ecological concerns and the current deforestation pressure faced by these ecosystems make them relevant biodiversity hotspots for conservation (Tse-ring et al. 2010).
The genus Eremanthus is composed of 22 tree and shrub species restricted to the arid Cerrado in the Central Plateau of Brazil (Loeuille et al. 2012). E. erythropappus (DC.) MacLeish occurs generally at altitudes of 700-2400 m as colonies within secondary forest, mainly in seasonally dry tropical forest of the coastal range and it is also found through the Brazilian savanna and the Atlantic forest. In these sites, the species has a key role in environmental services mainly as a source of water. Meanwhile, E. incanus (Less.) Less is particularly common in Minas Gerais State, at elevations from 800 to 1850 m in savanna, secondary forest, and semi-arid region. Possible hybrids of these two species have been found in southern Minas Gerais (MacLeish 1987).
The exploitation of these species dates from the Brazilian colonial period. At that time, its wood began to be used as fences material due to the long-standing durability, caused by the high concentration of essential oils. Because of its composition, E. erythropappus has been used for the extraction of alfa-bisobolol, an alcohol distilled from the essential oil extracted from its wood that is in high demand by pharmaceutical industries. However, the alfa-bisobolol extracted from the wood of E. incanus has lower quantity and quality; therefore, its main commercial use is the production of logs, largely for fence construction (Scolforo et al. 2012).
Several studies have analyzed the phenology and reproductive biology , genetic patterns (Barreira et al. 2006;Silva et al. 2007;Pádua et al. 2016) and in vitro propagation (Miranda et al. 2016) of Eremanthus species. The studies of its phylogeographic patterns and climatic variation responses are important for understanding historical patterns of seed dispersal, recolonization routes, genetic diversity distribution and species refugia.
Understanding the association between phylogeographical modeling and molecular data is an important method to better comprehend species' spatial dynamics (Collevatti et al. 2012). Besides, phylogeography examines the relationships between genetics and geographical species distribution and it is based on the spatial distribution of gene genealogies (Avise et al. 1987). The variability of conserved cytoplasmic genomes with uniparental inheritance, such as chloroplast for plants (cpDNA), results in low levels of mutation and no recombination making this marker suitable for phylogeographic studies (Avise 1994;Birky and Walsh 1988).
In this study, we evaluated the effects of climate change on the potential geographical distribution of E. erythropappus and E. incanus. We used chloroplast-specific markers to access the level and distribution of genetic diversity of E. erythropappus and E. incanus throughout Minas Gerais State, Brazil. We also used ecological niche modelling (ENM) to assess past and present distribution of these species and evaluate the possible impacts of future climate change on the species distribution for the year 2070. Our aim is to infer the spatial-temporal dynamics of these species, as suggested by predictions of past, current, and future species distribution.  Table 1. Sampled individuals were at least 50 m apart and each individual was georeferenced. Voucher specimens were deposited in the ESAL herbarium of the Federal University of Lavras (UFLA), Brazil.

RFLP analysis
The PCR products (5 lL) were digested for 1 h using the HinfI, HindIII, HaeIII, MseI, TaqI, AluI, EcoRI, and BamHI restriction enzymes. Each restriction was undertaken in a total volume of 10 lL at the temperatures recommended by the enzyme manufacturer (Gibco-BRL, Brazil). The digestion products were separated on 3% agarose gels for 2 h at 60 V and the gels were stained with ethidium bromide and photographed under UV light. The replicability for each PCR and enzyme/PCR product combination was assessed by replicating each reaction.

Genetic analysis
The total genetic diversity (H T and V T ), intrapopulation genetic diversity (H S and V S ), and the differentiation indexes (G ST and N ST ) were calculated based on alleles frequencies and the distance between haplotypes for nonordered and ordered alleles, respectively, using the software HAPLONST (Pons and Petit 1996). The same program was used to estimate v-type parameters. A minimum spanning tree (Excoffier and Smouse 1994) and AMOVA were calculated using ARLEQUIN 3.11 (Excoffier et al. 2007). The contribution of each population to total diversity and total allelic richness (CT) was estimated using the software CONTRIB (Petit et al. 1998). This contribution is explained by two factors: the diversity of the population (CS); and the differentiation of the remaining populations (CD).

Species distribution modeling
The database consisted of a set of 19 raster files of bioclimatic variables and altitude were extracted from WorldClim (www.worldclim.org/bioclim) (Hijmans et al. 2005). We used 341 herbarium occurrence locations of E. erythropappus (198) and E. incanus (143), downloaded from the Global Biodiversity Information Facility (GBIF 2018). The modeling was performed using the maximum entropy model software, MaxEnt 3.4 (Phillips 2017). According to Hernandez et al. (2006), this model has a high level of accuracy for species modeling and it is less sensitive to small sample sizes (Wisz et al. 2008). We used a spatial resolution of 2.5 arc-seconds (5 km at the equator). The database was divided in 25% for validation and 75% for training and to train the five periods under study: Last Glacial Maximum (LGM), about 22 kyr BP; Mid-Holocene, about 6 kyr BP; current conditions, between 1960 and 1990; RCP 2.6 and RCP 8.5 for the year 2070. RCP 2.6 represents the most optimistic scenario, with global temperature increase varying from 0.3 to 1.7°C by the year 2100, and RCP 8.5 represents the most pessimistic scenario, with an increase in temperature varying from 1.4 to 4.8°C by 2100 (Stocker et al. 2013).
We tested the Community Climate System Model (CCSM 4.0) (CESM 2017) to simulate the earth's climate conditions. This model has been widely used to understand past conditions and project the effects of future climate change (Gent 2011). We used 20 replicate cross validation runs, in which all occurrence data were randomly split into a number of equal-sized groups (Phillips and Dudík 2008). In addition, we used a maximum of 2000 iterations and convergence threshold of 0.001 (Vieira et al. 2015). The model efficiency was tested using the area under the ROC curve (AUC) generated by the MaxEnt model. The model values vary from 0 to 1, with values close to 1 showing a high efficiency and it is often used on validation of models (Fielding and Bell 1997;Elith et al. 2006;Phillips 2017).

Genetic analyses
An initial screening of 15 pairs of universal chloroplast primers with eight restriction enzymes revealed easily scorable variation for only one primer-enzyme combination (SFM/HinfI) which was used throughout the study. Altogether, this combination revealed one restriction site mutation and two insertion/deletion mutations, constituting three haplotypes (Fig. S1). The minimum spanning tree showed that haplotype 1 is the oldest and diverges from six mutations and haplotypes 2 and 3 are divergent from two mutations (Fig. 1). Haplotype 1 appears in very low frequency in both species (1.6% in E. erythropappus and 2.5% in E. incanus). We found that haplotype 3 was the most common (95.2%) in E. erythropappus and haplotype 2 was the most common for E. incanus (90%). We did not find haplotype 2 in E. erythropappus populations located in the center-south region, but it appears with a frequency of 11% in the northern populations (Fig. S2).
Populations RP, JO, and FS contributed most to the total diversity (CS) in E. erythropappus (37%, 13% and 11%), and population FS (42%) contributed most for E. incanus due to the presence of haplotypes 2 and 3. The contribution of differentiation (CD) in all populations was very low (2% in MP, BA, CA and AI populations for E. erythropappus and 8% in ST population for E. incanus) which is probably due to the restricted number of haplotypes found for these species (Fig. S3). The predominance of only one haplotype in all populations of both species generated low estimated diversities for both unordered and ordered alleles (Table 2). No genetic substructure was found for these species and regions (G ST = 0.611 and N ST = 0.541), indicating that the identified genetic differentiation is not due to the distribution of haplotypes (Table S1).
For E. erythropappus, the total genetic diversity (H T and V T ) were both equal to 0.167, and the intrapopulational genetic diversity (H S and V S ) were 0.174 and 0.117, respectively. In general, we found a higher variability (V S = 0.350 and V T = 0.312) on the northern populations when compared with the southern ones (V S and V T-= 0,027). For E. incanus, H T and V T were equal to 0.210 and 0.209, and H S and V S equal to 0.220 and 0.217, respectively (Table S1). We confirmed this result with AMOVA, which shows that most of the genetic variation occurs between species (Table 2). Despite the low variation between populations within species, the genetic differentiation was significant. The results from AMOVA showed that populations of E. erythropappus are structured regionally and most of the genetic variation occurs within populations (86%). This result was confirmed by the Mantel test made between the pairs of populations which indicated a positive and significant correlation between distance geography and genetic distance (r = 0.33, p = 0.04). We found no genetic variation among regions for E. incanus, demonstrating a high genetic variation within populations (97.7%). This result was confirmed by the Mantel test, which indicated negative but not significant correlation between populations (r = -0.36, p = 0.91).

Species distribution modeling
From the initial 20 variables extracted from WorldClim, 10 were chosen as they had the greatest influence on the MaxEnt model (Table 3) The ENM performed adequately for all five study periods (Last Glacial Maximum; Mid-Holocene; current day; RCP 2.6; and RCP 8.5). The environmental variables that most influenced the MaxEnt model were altitude, mean temperature of coldest quarter and precipitation seasonality (Table 3). We observed a gradual reduction in suitable area for species distribution over time (Figs. 3c, d, 4c and d).
The main difference was from the LGM climate scenario, the oldest period, where the species distribution was Fig. 1 Minimum spanning network for the three haplotypes obtained from Eremanthus erythropappus and E. incanus. The observed traits refer to the number of mutations observed among haplotypes. The numbers 1, 2, and 3 refer to the haplotypes significantly wider in comparison to present day with a remarkable reduction in suitable area for both future scenarios.
Based on the ENM for current conditions, we found that the areas of greatest suitability for both species are located in the mid-west, northeast, and southeast regions of Brazil, mainly in Minas Gerais State (Fig. 2). According to Camolesi (2007), the area of E. erythropappus occurrence covers 457 municipalities in Minas Gerais, which corresponds to 34.38% of the State's area. This species also occurs throughout the southeastern part of the Brazilian Central Plateau, Brazilian coast, Cerrado and rupestrian fields in the Federal District, and the States of Goiás, Espírito Santo, São Paulo, and Rio de Janeiro, Brazil (Macleish 1987; Loeuille et al. 2012). On the other hand, while E. incanus occurs naturally in montane grasslands in Bahia, Espírito Santo, Rio de Janeiro, and São Paulo, the species occurs abundantly in Minas Gerais State, forming large aggregates in the savanna and highlands (Araújo 1944).

Discussion
Haplotype Quaternary climatic fluctuations may have played an important role in shaping haplotype distribution for both E. erythropappus and E. incanus, which show similar patterns of phylogeographic distribution and respond similarly to the effects of climate change over time. Based on our analyses, the most intense haplotype diversification and gene flow probably occurred during the last glacial period when both species were distributed over a vast area in South America. This wide dispersal may be responsible for the occasional presence of the same haplotype in both species. The ENM predicted a reduction in suitable area during the Mid-Holocene and consequently a decrease in gene flow between populations for E. erythropappus and E. incanus. Additionally, ENM predicted an even greater reduction in suitable areas for the species for both studied future scenarios.

Genetic analyses
E. erythropappus and E. incanus populations were characterized by low cpDNA variation and overall diversity and no population differentiation. The results for E. erythropappus showed contrasting patterns to those reported in previous studies using nuclear markers, such as isozymes, RAPD and ISSR; these studies showed high population diversity and overall diversity and significant differentiation between populations (Estopa et al. 2006;Moura 2005;Pádua et al. 2016) probably due the higher mutation rate of those markers than cpDNA (Provan et al. 1999;Rawat et al. 2014). Unlike the cytoplasmic genome, the chloroplast genome is clonal, or inherited maternally, and historical events, such as glaciation and climatic change, have a greater impact on its variation over geological time (Avise 1994;Avise et al. 1987). The absence of population differentiation based on cpDNA markers for these species is probably due to their wide distribution during the LGM. Meanwhile, recent fragmentation and exploitation has led to population divergence as observed in studies using other molecular markers. The divergence between haplotypes 2 and 3 could have occurred alongside the separation of E. erythropappus and E. incanus into different species or they could represent older polymorphisms, with each haplotype being initially only partly associated with the species. These issues can be clarified with an increase in available cpDNA markers (e.g., by sequencing non-coding regions of the chloroplast genome) and by analyzing a larger sample of populations, species, and locations. With such data, we can better understand the taxonomic and geographical distribution of cpDNA variation for these species and their past evolutionary dynamics.
Haplotype 1 was found in only one E. erythropappus individual from the LV and AD populations, located in the southern and northern regions of Minas Gerais, respectively. Moreover, haplotype 2 was found in only three of the four E. erythropappus populations in the northern region of the state (PR, FS, and JO). Similarly, haplotype 1 was rare for the E. incanus populations, and only identified in one individual from the population of São Tomé das Letras (ST). In general, we observed a lower number of haplotypes for the populations from the central-south (EECS) and southern (EIS) regions for both species. Cloutier et al. (2005) studying 11 different populations of Carapa guianensis using chloroplast DNA variation, found low intrapopulational diversity for (H S and V S-= 0.02) and high values of total diversity (H T = 0.82 and V T = 0.73), and a population differentiation value (N ST ) equal to 0.97. On the other hand, Schlögl et al. (2007) studying the genetic diversity of Araucaria angustifolia, also using cpDNA markers, found values of H S and H T

Fig. 2 Geographical distribution of haplotypes and current distribution for (a) Eremanthus erythropappus and (b) Eremanthus incanus.
Colors indicate different haplotypes according to the legend. The circle size represents the sample size in each population and circle sections represent the haplotype frequency in each sampled population equal to 0.441 and 0.612 respectively, and a low N ST value (0.280). In our study, the low percentage of total diversity observed for E. erythropappus and E. incanus may be related to the small number of haplotypes. The populations with distinct haplotypes should be used as reference areas for ex situ conservation when moving individuals between areas or transplanting them to a new site (Chen et al. 2015).

Ecological niche modeling
The current distribution and genetic diversity of plants are the result of climatic events, such as the Pleistocene period of glaciation (Abbott and Brochmann 2003;Comes and Kadereit 1998;Hewitt 2000;Newton et al. 1999). The changes that occurred during the glacial period contributed to the contemporary spatial distribution of species through extinction, colonization, dispersals, and isolation (Abbott and Comes 2004;Newton et al. 1999;Lexer and Widmer 2001).
During the last glacial period, the Neotropical environment was both cooler and drier than current climatic conditions. In this environment, the vegetation, such as the seasonally dry tropical forest and the Cerrado, would have reflected those climatic conditions (Pennington et al. 2000). Prado and Gibbs (1993) concluded that fragmentary and mostly disconnected species distribution patterns of dry tropical forest are remnants of a once extensive and largely contiguous seasonal woodland formation that may have reached its maximum extension during a dry-cool period between 18 and 12 kyr BP, coinciding with the contraction of the humid forest.
Paleontological evidence from Á guas Emendadas, Brazil, dated to the period between 26 and 21.5 kyr BP, shows that the Cerrado biome was mainly characterized by abundant vegetation of grass, savanna, and gallery forests (Barberi et al. 2000) and dominated by Poaceae, Borreria, Cyperaceae, and Asteraceae species (Behling and Hooghiemstra 2001). Surprisingly, the currently occurrence of E. erythropappus and E. incanus is essentially situated in open, high-elevation, and grassland areas and it is restricted to either elevation over 1000 m or higher latitudes (Mayle et al. 2009).
We found a similar distribution pattern for E. erythropappus and E. incanus for past and current climatic conditions, thus reinforcing these theories for the LGM period (* 21 kyr BP). We found an intense distribution for both species, mainly in the northeast and southeast of Brazil, which largely represents the current Cerrado and some parts of the Caatinga biome (Figs. 3a and 4a). Additionally, the altitudes of occurrence of the studied species are mainly around 1000 m (Lorenzi 1992). In fact, this pattern of wide distribution for both species in the past may be responsible for the sharing of haplotypes between individuals of the northern (EEN) and southern regions (EECS), as supported by our genetic analyses.
Between 10 and 20 kyr BP, the area occupied by tropical rainforest in South America retreated as a result of climate change and previously continuous populations suffered a reduction in size and likely experienced genetic isolation among distant populations (De-Granville 1988;Prance 1973Prance , 1982. For these species, the loss of genetic diversity within populations and the increase in spatial genetic structure are expected to produce a bottleneck effect (Nei et al. 1975). Our results for ENM may suggest a similar pattern of a significant reduction in species distribution areas for the Mid-Holocene (Figs. 3b and 4b) and for the current conditions ( Fig. 2a and b).
The vast area previously occupied by both species during the LGM (Figs. 3a and 4a) appears have offered the appropriate environmental characteristics for the species to share haplotypes and may have served as potential refugia for tropical trees. This is likely the result of greater connectivity between populations and, consequently, more favorable conditions for the species to exchange alleles. Thus, the reduction in the suitable area over time may be responsible for the separation of these populations, although currently maintaining few haplotypes. A similar pattern as found by Carvalho et al. (2017) with an ENM analysis of the distribution and genetic diversity of Euterpe edulis species in the Brazilian Atlantic coast, indicating a reduction of climatically suitable areas over the time.
According to Behling and Hooghiemstra (2001), the Cerrado domain only became fully established after 7 kyr BP. However, throughout the early Holocene, from 9.7 to 5.5 kyr BP, the Cerrado was already reflecting a pattern of low precipitation and a dry season of 6 months, with the formation of some gallery forest. Furthermore, between 6 and 5 kyr BP, the climate in most South American savannas was drier than during the late glacial and late Holocene periods, while distribution of the savanna (8.2 kyr BP) in the early Holocene was much wider than during the late Holocene (4.2 kyr BP) (Behling 2001). Similarly, our findings suggest that both species had experienced a reduction beginning in the last glacial period until current conditions, probably due to substantial forest reduction and expansion of the savannas during the mid-Holocene (5 kyr BP) (Ledru 1993).
Because it occurs on fertile soils, seasonally dry tropical forest areas have been widely converted to agricultural land use, which has resulted in the vast destruction of these forests in many Neotropical dry forest biomes (Murphy and Lugo 1995;Pennington et al. 2000). Additionally, global climate change may have a severe impact on the functioning, distribution, and establishment of species in ecosystems (Nabout et al. 2016). ENM has been widely used to predict the effects of climate change on species Effects of climate changes on distribution of Eremanthus erythropappus and E. incanus … 359 distribution (Beck 2013;Peterson 2011). ENM for current and future projections showed a higher probability of occurrence for both species in the southern region of Minas Gerais State. This pattern can be the result of a reduction in suitable areas in the northern part of the state over time.
The current potential geographic distribution of E. erythropappus and E. incanus shows a wide distribution in the Cerrado and Caatinga biomes. However, our results indicate that the future climate change scenarios studied herein will result in a significant reduction in the distribution of both studied species (Figs. 3c,d,4c,d). Our study highlights the potential of bioclimatic variables to predict the effects of climate change for conservation purposes. We recommend the development of public policies that aim to include species genetic attributes in the maintenance of species viability. Such policies are necessary to ensure that both species studied herein will continue to reproduce and maintain diversity for future generations. Furthermore, the

Conclusion
Priority areas for the conservation of E. erythropapus include the Medina population due to the presence of haplotype 1 and Franscisco Sá, Joaíma, and Rio Pardo de Minas due to the occurrence of haplotype 2, which is exclusive to the northern region of Minas Gerais State. In the center-south region, the population of Lavras is also a priority because of the presence of haplotype 1, which differs from the other populations in the region. For E. incanus, it is necessary to preserve the populations of Franscisco Sá (north) and Lavras (center-south) as they present haplotype 3, and São Tomé das Letras because it is the only population to show haplotype 1. The ecological niche modeling revealed an extensive area of occurrence of E. erythropappus and E. incanus during the last glacial maximum and also indicated a reduction in distribution following both the Mid-Holocene and current conditions. Additionally, the predictions for the year 2070 showed an even greater reduction in suitable area as a result of climate change.