Future uncertainties for the distribution and conservation of Paubrasilia echinata under climate change

Paubrasilia echinata is a widely cultivated endangered tree species with small populations restricted to a narrow strip of habitats along the Brazilian coast. The potential impacts of climate change on the distribution of P. echinata have yet to be investigated, and so it remains unknown whether protected areas will ensure the persistence of the species in the future. Here, we estimate the impacts of climate change on the distribution of P. echinata inside and outside of protected areas considering different climate change scenarios and two different sets of presence records: natural distribution and cultivated records. Future scenarios showed a gradual reduction in climatically suitable area both inside and outside of protected areas. Projections indicate a trend for a shift to the highlands of Southeast Brazil, and the loss of several areas throughout the entire distribution of the species. Predicted climatic conditions will be unsuitable for P. echinata inside most protected areas. Information provided here will be relevant in planning future national actions for this species, which is a must to properly protect this long-exploited tree.


Introduction
Climate change is a major threat to biodiversity and natural ecosystems on a global scale (Parmesan & Yohe 2003;Root et al. 2003). Shifts in species phenology and distribution in response to ongoing climate change have been reported for marine, freshwater, and terrestrial groups (Parmesan 2006;Franklin et al. 2016). Several species of vertebrates and plants are expected to lose suitable climate conditions both inside and outside of protected areas (Araújo et al. 2011;Li et al. 2015). Conserving future biodiversity in a changing climate is a complex issue, as some protected areas will become less effective in the future, while others will maintain their conservation value (Thomas & Gillingham 2015;Berteaux et al. 2018). Therefore, predicting future climate scenarios and their effects on the species ranges is crucial to effective conservation strategies and to mitigate the impacts of climate change on biodiversity.
When planning for conservation strategies, it is crucial to take climate change and species redistribution into account, if not, jeopardizing species persistence. Ecological niche modeling is the main tool to predict the impacts of climate change in species distribution (Peterson et al. 2000;2002;Peterson 2001). With the development of new algorithms (Guisan & Thuiller 2005) and more accurate future climate predictions with multiple scenarios (Moss et al. 2010), we can reduce uncertainties to better support conservation decisions (Elith et al. 2006).
Paubrasilia echinata (Fabaceae), locally known as 'Pau-Brasil' is an endemic and endangered tree species with small populations restricted to narrow habitats along the Brazilian coast. Pau-Brasil is estimated to have been widely distributed before the arrival of Europeans, as documented by Américo Vespúcio (Fontana et al. 1994) when he wrote to the king of Portugal saying the Brazilian coast had no profitable resource, but an infinity of Pau-Brasil. Commercial exploitation until 1875 drastically reduced populations (Souza 1939), despite conservation attempts such as the Pau-Brasil statute (Brasil 1605). Since the development of synthetic dyes, the main use of Pau-Brasil shifted to luthiery, as the building material of highend violin and cello bows (Skeaping 1955;Longui et al. 2010). In 1978, Pau-Brasil was declared as Brazil's national tree (Lei N° 6.607). Currently it is widely distributed due to conservation efforts. It is cultivated in plant nurseries for urban landscaping and educational purposes. Although the species is found in 11 protected areas, it is unknown whether these areas will maintain viable populations under global climate change, which is the main driver of species redistribution in the XXI century (Thuiller 2004;Diniz-Filho et al. 2009;Thuiller et al. 2011;Loyola et al. 2014;Sales et al. 2017). So far no studies have estimated the potential impacts of climate change on the distribution of P. echinata. This is particularly alarming given the low levels of gene flow and genetic diversity of its populations in nature (Lira et al. 2003;Cardoso et al. 2005).
Our goal is to estimate the impacts of climate changes on the distribution of P. echinata, searching for climatic stable areas where a long-term conservation strategy could be applied, as well as to assess its potential distribution inside of protected areas in different climate change scenarios. We expect a decline of climatically suitable areas for P. echinata both inside and outside of protected areas, as the Atlantic forest habitats where it occurs are highly vulnerable to climate changes (Bellard et al. 2014;Esser et al. 2019). A secondary goal is to estimate the contribution of cultivated records to the models. We expect those records to help find sites where the species will have a long-term persistence and in which the species would never reach despite human mediated migration. It is important to point out that the use of non-natural records is not recommended in ecological niche modeling, since they may include biased information (i.e. environmental information where the species is not actually suitable without human management), and further inferences over this results should be taken with caution.

Study species
Paubrasilia echinata (Lam.) Gagnon, H.C.Lima & G.P.Lewis is a semideciduous or deciduous medium sized tree (up to 20m tall), armed with prickles, chestnut brown to almost black bark; red heartwood and when injured the trunk exudates a red sap (Gagnon et al. 2016). Paubrasilia is a monospecific genus endemic to Eastern Brazil, with occurrence from the Brazilian states of Rio Grande do Norte to Rio de Janeiro within the Brazilian Atlantic Forest hotspot. Its main habitats are well-drained soils in coastal white-sand woodlands (tall restinga), seasonal forest (semideciduous forest) and rain forest (Gagnon et al. 2016;Lima 2019). This species is listed in CITES Appendix II and is listed as endangered at the national (CNCFlora 2012) and global (IUCN 2019) level due to a population size reduction of 50 % in three generations.

Datasets
Species occurrence data was obtained from SpeciesLink, GBIF, JABOT and Reflora, and were subject to a data cleaning approach, where we deleted miswritten coordinates, as well as records with zero degrees. We also deleted duplicated records within the same cell of a raster with 5-arcmin resolution. Data were then divided in a natural distribution and a cultivated group. The first group was build considering Castro (2002), which after an extensive literature review, based on 23 studies, proposed that the natural distribution of P. echinata is the brazilian coast, extending from 5º39'S in Rio Grande do Norte to 23 o S in Rio de Janeiro. Records selection was made manually using QGIS. The second group comprised all records, including both natural distribution and cultivated records (Fig. 1). Presence records summed a total of 237 sites, of which 180 were considered from its natural distribution (Tabs. S1, S2 in supplementary material). We considered this second approach due to the potential resistance to be unveiled by environmental conditions in those sites, expanding the niche breath of the species. The following routine was then applied for both datasets.

Variable selection
Variable selection was conducted in two steps. We firstly pre-selected WorldClim 1.4 (Hijmans et al. 2006) 5-arcmin bioclimatic variables that matched species ecology and habitats where it is found (Esser et al. 2019). Rainforest variables were mean temperature of wettest quarter (BIO 8) and precipitation of warmest quarter (BIO 18), which are both proxy of tropical climate, and annual precipitation (BIO 12), which, together with BIO 18, tracks the lack of dry season. Restinga forest was represented by mean diurnal range, temperature seasonality and precipitation of driest month (BIO 2, 4 and 14), which are respectively proxy of maritimity, altitude and dry. BIO 14 was also representative of seasonal forest, since it is also a dry formation. To represent this last habitat, were selected also mean temperature of driest quarter (BIO 9), which is proxy of dry and heat, and precipitation seasonality (BIO 15), which is proxy of seasonality. The summarized pre-selection is presented in Table 1. Secondly, we conducted a statistical selection using a Variance Inflation Factor (VIF) approach with usdm package from R environment (Naimi et al. 2014; R Development Core Team 2018), selecting variables with correlation between 0.5 and -0.5; and VIF lower than 3. This grants a lower dependency from variables and high predictability. This statistical approach is dependent on the background data, so we masked bioclimatic variables with a one-degree width buffer for each presence coordinate, which better describes the distribution of the species and their surroundings. Variables marked with a * are those that most appeared in the selection routine, thus used to build models. . Finally, we calculated the committee average of projections, which returns both a prediction and a measure of uncertainty, since it is the simple average of binary predictions (i.e. cells with values close to one are those that most models agree with a presence, cells with values close to zero are those that most models agree with an absence and cells with value of 0.5 are those that half of models predict a presence, while the other half predict an absence). Ecological niche models were made using biomod2 package in R.

Area calculation
Area value was obtained with a weighted method, where presence probability is multiplied by cell's area, then summing all raster values (Esser et al. 2019). This approach results in a conservative area that considers that occupancy is not equal between cells. When we first obtained area values we realized the Amazônia biome had a great bias in results. Given the fact that its soil properties are a key factor to species persistence and we have not considered it in these models, we decided to exclude suitability values in Amazônia in this step. Despite this fact, we present here the complete projections.

Results
Models TSS' values ranged in natural distribution from 0.186 to 0.744 with mean value of 0.562 and standard deviation of 0.112. When considering complete data, the values ranged from 0.151 to 0.787 with mean value of 0.569 and standard deviation of 0.116. Ecological niche model unveiled that currently P. echinata has a potential suitable distribution ranging from 1,137,920 to 1,794,434 km 2 (Tab. 2). Of this total area, 7.10 % to 7.39 % are considered to be within protected areas (Tab. 3). Future scenarios showed a gradual reduction in species suitable area (with increasing CO 2 concentration), ranging from 1,648,148 km 2 for RCP 2.6 in the year 2070, considering both natural and cultivated records, to 376,479 km 2 for RCP 8.5 in the year 2070, considering just records in the natural occurrence area (Tab. 2). In the same way, potential distribution inside protected areas in future scenarios ranges from 144,155 km 2 for RCP 4.5 in the year 2070, considering both natural and cultivated records, to 45,645 km 2 for RCP 8.5 in the year 2070, considering records in the natural occurrence area (Tab. 3).  Projections from natural occurrence models revealed a considerably larger distribution than the previously known for P. echinata (Fig. 2). The states of Mato Grosso, Goiás, Minas Gerais and the continental portion of Bahia seem to have the same environmental suitability as its natural distribution. Furthermore, the state of São Paulo, as well as scattered sites in Mato Grosso do Sul and even Paraná may also be considered good places for that species to grow.
Future projections indicate a shift trend to Southeast Brazil and loss of several areas throughout its distribution range. It is expected to remain, in the worst cases, in few places along the Brazilian coast (Espírito Santo and Alagoas) and in the state of Minas Gerais and São Paulo (Fig. 3). Considering cultivated records, projections showed a similar trend along its distribution, but with higher overall suitability, remaining, in the worst case, with suitable areas in the states of São Paulo and Rio de Janeiro in addition to projections considering just natural occurrences (Fig. 4). Most of all, the coastal region of Pernambuco and Alagoas, as well as the South of Espírito Santo and São Paulo are key regions to the preservation of P. echinata, since they will have high suitability even in the worst scenarios. We identified three main climatic stable areas: the Northeast Brazil (Fig.  5A, a), the highlands of Southeast Brazil (mountains of Espírito Santo and Minas Gerais, Fig. 5A, b) and highlands of São Paulo (Fig. 5A, c). Stable areas presented a low number of future scenarios agreeing with the presence of P. echinata inside most protected areas.

Discussion
Climatically suitable area for P. echinata is larger than previously considered (Gagnon et al. 2016), although it may not fully colonize area within the continent since the coast of Brazil has a major altitudinal shift. In this way, sierras of the Brazilian coast are possibly acting as a barrier to migration, even though there is a suitable area for the species to grow. The Mean Diurnal Temperature Range (BIO2) variable was firstly selected to generate models, since it is a proxy of continentality, but it had high correlation with other variables, thus it was discarded after our selection routine. P. echinata thrives in a relatively wide variety of habitats, from coastal outcrops and white-sand forests to seasonal and rain forests, but always in well-drained soils (Gagnon et al. 2016). To consider soil attributes in ecological niche models may be particularly difficult, especially when projecting models into future scenarios, since soil will not change as quickly as the climate, resulting in a higher importance to this variable. The lack of studies on the species distribution limits our inferences over which environmental filters are working to prevent P. echinata's expansion into areas where it could potentially occur. Its relation with well-drained soils may give us a clue on why it does not colonize inner continent areas and other biomes such as Amazônia and Pantanal, despite the evidenced climatic suitability. Regardless of that, projections show an altitudinal shift in suitability through future scenarios.
As we didn't consider landscape metrics in this study, such as structural and functional connectivity, edge effect or patch size, our inferences over populations dynamics are weak, and it is essential in further studies. The increase of relatively protected area within future scenarios, followed by a decrease in absolute protected and total area, indicates that habitat loss will be more intense outside conservation units. This leads us to a problematic scenario, since a major  driver of species extinction is habitat fragmentation, which is particularly severe in P. echinata populations due to exploitation and selective logging (Cardoso et al. 2005). Despite that, it makes necessary to point out that, due to this long-term relative isolation of the natural populations, there is a high genetic diversity, with markedly morphological differences that may be recognized in future studies as new taxa (Lima 2019), and which may grant species quick adaptation, resilience and persistence through scenarios, although genetic diversity within populations is very low due to the discontinuity of the distribution.
Despite the fact that none of the protected areas that houses Pau-brasil nowadays may sustain it in future scenarios, the species will not be unprovided of protection. Parque Nacional da Serra de Itabaiana (Sergipe), Refúgio da Vida Silvestre Mata do Urucu (Pernambuco), Reserva Biológica da Pedra Talhada (Alagoas), Reserva Biológica Augusto Ruschi (Espírito Santo), Área de Proteção Ambiental Bacia do Paraíba do Sul (São Paulo) and Parque Estadual da Serra do Mar (São Paulo), as well as three different private reserves (Reservas Particulares do Patrimônio Natural -RPPNs) are predicted to harbor P. echinata in the next 50 years.
Climate stability maps (Fig. 5) can also be interpreted as potential areas for reforestation and restoration, which is particularly interesting for P. echinata, since light intensity in open areas of the forest positively influences growth rates (Zani et al. 2012). This species can act as a colonizer and thus as a promoter of the restoration process in degraded sections on forest edges (Sattler et al. 2018), facilitating as well the species migration.
Paubrasilia echinata will have a major obstacle in future scenarios, which is to adapt to new environments where it is not growing naturally. It may be capable to do that, since cultivated organisms show a high adaptation capacity to a myriad of habitats, developing relatively well. The significant difference in area between models (natural vs. cultivated) unveiled a huge potential for species adaptation. To consider cultivated records can drive future researches with ecological niche models to other cultivated species as well. Although some cultivated organisms may not bloom or fruit (e.g. due to ecological needs) or generate viable seeds (e.g. due to pollinators absence), these records amplify the niche breath of the species.
Even if P. echinata suffers greatly with climate changes, it is not likely to be extinct since it is an important symbol of our country, having not just a crucial role in Brazil's history, but also in the development of Brazilian culture. P. echinata may be a potential flagship species for the conservation of Brazilian flora, since it has public appeal and is not restricted to one habitat. This research may also be key to a future national action plan for this species, which is a must to properly protect this long-exploited tree.