Due to the highfrequency of severe weather events, changing climate has recently become matter of concern across the globe. Global mean surface temperatures have increased 1.07°C till 2019 since the industrial revolution and expected to increase by 3.3°C to 5.7°C under extreme scenario (shared socio-economic pathways [SSP]5-8.5) by the end of century (Intergovernmental Panel on Climate Change [IPCC] 2021). Nearly one-quarter of all species of plants are thought to be endangered as per the IPCC’s fifth assessment report (Stocker 2013). Climate change, which is currently altering species’ ranges and will continue to do so in future at more severe rates, will render harm to all restoration efforts (Weiskopf et al. 2020). Restoration strategies need to integrate novel scientific techniques like species distribution models (SDMs) and Nature Based Solutions to fulfill the requirements of current and future landscape issues (Beatty et al. 2018).

SDMs are being used to forecast the future distribution shifts caused by changes in climatic factors or land-use changes (Guisan and Thuiller 2005). Conservationists have used SDMs to anticipate future landscapes by employing factors that are either directly or indirectly linked to recently accessed global climate data (Pereira et al. 2010). SDMs can help in identifying corridors between conservation areas which can facilitate mobility across temperature differences, and can help in making strategy for enhanced extreme weather events in attempt to conserve a threatened species (Bateman et al. 2012; Nuñez et al. 2013). SDMs are quite useful tool to plan afforestation projects by projecting suitable areas for species under different climate change pathways (Hidalgo et al. 2008). Afforestation projects even in degraded landscapes of the world have been reported as successful, however, climatic factors such as climate extremes, seasonality factors and species selection play a vital role in success of any afforestation project (Cuong et al. 2019; Li et al. 2016; Majumder et al. 2013). High temperatures, drought stress, disease and low soil nutrient levels are factors responsible for extensive dieback and mortality of many tree species (Allen et al. 2010; Colangelo et al. 2018; Gazol and Camarero 2022).

Many afforestation projects are being carried out all around the globe, without considering niche dynamics of planted species. Due to this negligence of concern governments, a lot of species were introduced in their non-native zones (Li et al. 2018). The climate has been found as a significant factor influencing the large-scale distribution of many species (Poortinga et al. 2019). Over the last 30 years, global climate change has been shown to cause alterations in the distribution of many species and it may be the dominating driver contributing ultimately to species extinction in the short term or it may have a synergistic implication with other extinction causes (Román-Palacios and Wiens 2020). Many studies have predicted potential shifts in the tree line, changes in distribution of tree species, variations in forest type borders and interspecific competition under different climate scenarios (Boisvert-Marsh et al. 2014, 2019; Iverson et al. 2019; Meng et al. 2021; Zhu et al. 2018). Taking these into account, it is evident that climatic factors and species selection are important factors for afforestation projects. According to the Bonn challenge, the global effort is to restore the 150 million hectares of deforested land by 2020 and further 350 hectares till 2030 launched by Germany and the International Union for Conservation of Nature (IUCN) in Bonn, 2011. By meeting these targets, yearly emissions can be cut up to 8.8 billion tons by 2030 (UN Climate Summit 2014). The Bonn Challenge amplifies current goals of Global Partnerships on Forest Landscape Restoration; however, it also supports Aichi targets of the Convention on Biological Diversity (CBD) and for Reduce Emissions from Deforestation and forest Degradation (REDD+) (Pistorius and Freiberg 2014; Pistorius et al. 2017a, 2017b).

In view of Bonn Challenge, Billion Tree Tsunami Afforestation Project (BTAP) was initiated by Khyber Pakhtunkhwa (KP) Government in Pakistan, through KP Forest Department. Intensive restoration work was carried out in three forest regions with combination of 60% regeneration and 40% designed afforestation plantation. This afforestation project was completed in two phases (World Wide Fund for Nature Pakistan 2017). However, little is known about the criteria for site selection and species selection for the afforestation project (Sabir et al. 2022). Therefore, the goal of the present research is to determine the distribution characteristics of the habitat for the afforested species in present and future climate scenarios used for afforestation project BTAP in KP, Pakistan through ecological niche modelling. Ecological niche modelling of tree is beneficial to understand the relationship between species and their suitable habitat, predicting the responses to climate change and for sustainable afforestation programs (Carvalho et al. 2017). The specific objectives of this research were: (i) to identify significant change in the potential distribution of eight tree species of KP, Pakistan, used in afforestation project as a function of 19 climatic variables as well as soil, irrigation and elevation, (ii) to identify species with risks of decreasing their potential distribution area in the future climate under SSP2-4.5 and SSP5-8.5 for years 2041–2060 and 2081–2100, (iii) to verify the afforestation plantations sites of eight tree species, that is, if they fall within the predicted potential distributions under current and future climate.

Study area

The study area of this research is KP which is located in the north-west of Pakistan (Fig. 1). It covers 101,741 km² area, smallest province of Pakistan with respect to land. Wide altitudinal differences in the province result in various ecological climates ranging from 250 m a.s.l. in the southern part of province to 7,708 m a.s.l. in the northern part. The temperature varies from –14°C in the north to 51°C in the south. The annual precipitation varies from 130 mm in the south to 3,200 mm in the north (Ali and Begum 2015). KP contains Himalayan dry temperate forests between 1,525 to 3,350 m a.s.l. mostly present in inner Himalayas of North KP outside monsoon range, Cedrus deodara, Pinus wallichiana, and Pinus gerardiana as the main species. Himalayan moist temperate forests are located between 1,375 to 3,050 m a.s.l., which arethe outer Himalayas containing monsoon rainfall range. These contain species like, P. wallichiana, Abies webiana, and Picea smithiana. Sub-tropical chir pine forests can be traced in lower region of Himalayas at 900 to 1,700 m a.s.l. 63.6% of these forests occurred in KP. Scrub forests occur at 400 to 1,000 m a.s.l. on lower slopes of Himalayas in KP. Olea ferruginea, Acacia modesta are the main species of scrub forests. Tropical thorn forests occur in plain areas of KP at 400-m a.s.l. Acacia nilotica and Acacia Senegala are found in these forests. Oak forests in KP can be found from 1,200 to 1,500 m a.s.l. in moist and dry temperate zones (FAO 2020).

Figure 1. Map of the study area.

Occurrence data

Occurrence data for the respective species C. deodara, P. wallichiana, Dalbergia sissoo, Eucalyptus camaldulensis, Populus ciliata, Juglans regia, Senegalia modesta, and Vachellia nilotica (Table 1) were gathered from open database (Global Biodiversity Information Facility [GBIF] http://www.gbif.org). These data records generally represent only presences of species that can be used by modelling methods like Maxent (Dennis and Thomas 2000). Erroneous identifications, missing one or both coordinates, set coordinates to zero, swapped coordinates, positional mistakes, and duplicate records were among the most common errors. These data points were filtered, duplicate records were eliminated, occurrences without coordinates were removed, occurrence sites were validated by plotting, and species records were visually examined on a Google mapand maintained data in separate “.txt” files for each species (Sillero et al. 2021). Further, duplicated function of R was used for identification and removal of duplicated records. A detailed field surveys all over KP (study area) were conducted and selected only data points which originally exists on field. Data points for E. camaldulensis and V. nilotica were collected through detailed survey and only took block plantation data points in saline and water-logged soil areas of KP.

Table 1 . Number of occurrence points used for each species and mean AUC along with standard deviation values.

Species name	Number of occurrences	Mean AUC values	Mean standard deviation	Feature class	Regularization multiplier
Cedrus deodara	50	0.863	± 0.043	LQH	1.25
Dalbergia sissoo	15	0.866	± 0.043	LQH	1.50
Juglans regia	15	0.703	± 0.146	LQH	1.25
Pinus wallichiana	15	0.768	± 0.163	LQH	1.50
Eucalyptus camaldulensis	15	0.891	± 0.084	LQH	1.50
Senegaliamodesta	16	0.772	± 0.094	LQH	1.25
Populus ciliata	30	0.765	± 0.073	LQH	1.25
Vachellia nilotica	15	0.820	± 0.110	LQ	1.25

AUC: area under the curve; LQH: linear quadratic hinge; LQ: linear quadratic.

Billion tree afforestation project data

Data regarding KP afforestation project BTAP was acquired from BTAP Office, Peshawar, KP.

Bioclimatic data

Bioclimatic data which contains 19 variables was downloaded from WorldClim database version 2.1 from (www.worldclim.org) in raster format at 2.5 minutes spatial resolution. For the years 2021–2040, 2041–2060, 2061–2080, and 2081–2100, WorldClimate provides 9 global climate models (GCMs) and 4 SSPs. The Scenario Model Inter comparison Project (Scenario MIP) is a novel climate change scenario presented as part of the sixth phase of the Coupled Model Inter comparison Project (CMIP6). Scenario MIP integrates many common socioeconomic paths with recent anthropogenic emissions trends (SSP-RCP) (Eyring et al. 2016). SSP1-2.6 proposes a green development trajectory, with a multi-model mean warming of less than 1.5°C by 2100.SSP2-4.5 is a middle development route with a radiation forcing of around 4.5 W/m² in 2100.SSP3-7.0 is characterized by a high level of societal vulnerability as well as a reasonably high level of human radiation forcing. With a radiation forcing of 8.5 W/m² in 2100, SSP5-8.5 represents the extreme end of the range of future routes (Jones et al. 2016).

For this research, we adopted GCM, IPSL-CM6A-LR from CMIP6 because of its low-biasness and high-sensitivity (Boucher et al. 2020). Two shared socioeconomic pathways were selected (O’Neill et al. 2017) to evaluate the contrasting patterns of possible and distinct future climate conditions: the SSP2-4.5 (ambitious) as well as SSP5-8.5 (realistic). To evaluate the climate change impact by both the middle and end of the 21st century, two periods were utilized for SSPs, 2041–2060 and 2081–2100. Other variables such as elevation (Shuttle Radar Topography Mission digital elevation model, SRTM DEM) and soil types were also integrated into the model. Elevation (DEM) data was downloaded from (www.worldclim.org) at 2.5 minutes spatial resolution and irrigation (Irrigation Areas v.5) and soil type, were obtained from Food and Agriculture Organization of the United Nations (Fischer et al. 2008). To prevent multi collinearity, generated a collinearity matrix by using the raster.cor.matrix function of environmental niche models Tools package and removed variables with a correlation value (r) of 0.80 or higher (Fig. 2) (Holzmann et al. 2015; Kramer-Schadt et al. 2013).

Figure 2. Correlation test for variables.

Maxent model

The bioclimatic data raster files for current climate and future scenarios were clipped into KP extent. The text file of occurrence data, 13 non-colinear variables including bioclimatic variables, elevation (DEM), and soil type in ASCII format were considered for model input. Soil type was taken as categorical variable while DEM was taken as continuous variable. The maximum number of background points was set as 5,000 and 5-folds cross validation were used, splitting the data into 5 training (80%) and testing (20%) subsets that used each records once as training or testing. Each k-fold was run and averaged for area under the curve (AUC) and receiver operating characteristic (ROC). Furthermore, non-collinear variables, which contributed ≤ 0.5 percent to the Maxent model in dry run for each species were excluded in the final run (Table 2). The regularization multiplier and feature classes are the two main parameters in Maxent. The regularization multiplier defines how centered or precisely the output distribution is fitted (Morales et al. 2017). A lower number than the default of 1.0 would produce amuch-restricted output distribution which satisfies the provided presence data. A greater regularization multiplier would result in a forecast which is more spread out and least localized. For this research, initially in dry runs model regularization were checked with 8 regularization multipliers (1–3) at intervals of 0.25 to get the best optimized for regularization number for each species (Table 1). Eucalyptus camaldulensis, P. wallichiana were operated on 1.5 regularization multiplier and rest of the species on 1.25.

Table 2 . Percent contribution of variables.

Variables	Cedrus deodara	Dalbergia sissoo	Juglan regia	Pinus wallichiana	Eucalyptus clamedius	Senegalia modesta	Populus ciliata	Vachellia nilotica
Annual mean temperature	0.6	11.9	16.3	3.4	0	0	20.3	16.3
Temperature seasonality	0.1	1	0	3.7	3.6	2.3	8.7	1
Temperature annual range	35.5	0.1	0	32	0	0.1	0.1	0
Annual precipitation	0	0.1	0	13.7	0	0	2.7	0
Precipitation of wettest month	0.6	0	0	0.3	0	10	6.7	5.6
Precipitation of driest month	0.2	7	12.8	5.3	1.8	2.5	0	0.2
Precipitation seasonality	0	0	6.2	0	0	0	0.1	0
Precipitation of wettest quarter	0.3	0	0	0.2	4.9	1.8	0.1	0.4
Precipitation of driest quarter	0.2	0.8	10.9	0	2.1	1.2	0	0
Precipitation of warmest quarter	1.7	13	0	0	7.3	26.5	12.5	2.1
Precipitation of coldest quarter	0.4	0	0	0.5	2.1	0	0	1.7
Soil type	53.6	27.2	53.6	20.1	47.7	29.8	38.8	49.5
Elevation	6.8	38.9	0.2	20.8	30.6	25.9	9.9	23.3

The term “feature class” refers to a bigger set of transformations applied to the original covariates (Elith et al. 2011). This restricts the calculated probability distribution. The amount of species occurrence points determines the feature class selection. If there are few data pointsavailable, the algorithm confines the model to basic features by default. The model generally useslinear feature howeverwith at least 10 samples, the quadratic feature is employed; with at least 15 samples, the hinge is employed; and with at least 80 samples, the threshold and product are employed (Li et al. 2020; Phillips 2017). In this research, grid search approach was used to find the optimum feature classes, several numbers of dry runs were compared to select the best combination of feature class and regularization multiplier (Li et al. 2020). Two feature classes were used in this study linear quadratic and linear quadratic hinge (Table 2). All species distribution modelling was performed using Maxent (version 3.4.1) in the R package dismo (Hijmans et al. 2011).

Comparable maps for the current and the two future scenarios were produced after Maxent model was simulated. For the 2 future scenarios, we generated maps for 2 different timescales (2041–2061 and 2081–2100). The AUC values from the ROC were examined in order to assess the model’s efficiency. With all practicable parameters, AUC determines the model’s performance. When the AUC value is around 0.5 and 0.6, the model accuracy is considered failed, bad when it is within 0.6 and 0.7, fair when it is within 0.7 and 0.8, good when it is within 0.8 and 0.9, and outstanding when it is within 0.9 and 1 (Araújo et al. 2005). The possible habitat was divided into 4 levels using the Jenks’ natural breaks method: unsuitable environment (0.0–0.1), low suitable habitat (0.1–0.3), moderate suitable habitat (0.3–0.5), and high suitable habitat (0.5–1).

By using a likelihood threshold (0.8) that reflected optimum specificity at optimum sensitivity (Liu et al. 2005), we turned each uniform terrain used in the aggregated projection into a binary projection. The binary projection that resulted, assigned a value of 1 and 0 per each 1 km² pixel, indicating whether it would be suitable or non-suitable. Additionally, we utilized BTAP afforestation sites coordinates with respect to species and plotted those coordinates on binary projections to identify if the afforested sites fall within the suitable area in all considered scenarios. The Maxent model in R was used to simulate and analyze the results.

Model performance

Models for all the species performed well with the given set of training and test data according to their mean AUC training values (Table 1).

Tree species distribution

The model for C. deodara predicted most of its suitable environment in the north of KP at higher altitudes in current climatic conditions. The distribution of C. deodara can be seen confined in same regions under current and future climate projections (Figs. 3A and 4A). For future SSP2-4.5 (2041–2060) scenario, there was sudden decrease in suitable area however, some reduction in northern part was also observed. For future SSP2-4.5 (2081–2100) scenario more suitable areas emerged in southwest where there was an overall expansion (Figs. 3B, C and 4B, C). In future scenario of SSP5-8.5 (2041–2060), the distribution of C. deodara decreased with slight loss in the northeast. In future SSP5-8.5 (2081–2100) scenario, a slight increase in suitability was observed along with expansion in northern area (Figs. 3D, E and 4D, E).

Figure 3. Areas ranging from suitable to non-suitable are depicted as raw-value maps derived from Maxent exponential output. (A) Current. (B) SSP2-4.5 (2041–2060). (C) SSP2-4.5 (2081–2100). (D) SSP5-8.5 (2041–2060). (E) SSP5-8.5 (2081–2100). SSP: shared socio-economic pathway.

Figure 4. Binary maps for suitable and non-suitable areas with threshold value of 0.8 from Maxent output. (A) Current. (B) SSP2-4.5 (2041–2060). (C) SSP2-4.5 (2081–2100). (D) SSP5-8.5 (2041–2060). (E) SSP5-8.5 (2081–2100). Yellow dots are recent afforestation sites for respective species. SSP: shared socio-economic pathway.

The model for D. sissoo predicted most of its suitable environment in south-east of KP in current climate conditions (Fig. 3A). According to binary map (Fig. 4A), the model predicted very low suitability areas for D. sissoo. In SSP2-4.5 (2041–2060) scenario, there was overall expansion however, some increase in southern part was observed. In SSP2-4.5 (2081–2100) scenario, distribution in existing environment was observed (Figs. 3B, C and 4B, C). In SSP5-8.5 (2041–2060) scenario, similar pattern was observed with slight migration towards north. In SSP5-8.5 (2081–2100) scenario, distribution similar pattern was observed along with expansion towards north (Figs. 3D, E and 4D, E).

The model for J. regia incurrent climatic conditions predicted most of its suitable environment in north-west of KP, however few patches of suitable environment were also predicted in south (Figs. 3A and 4A). For future scenario of SSP2-4.5 (2041–2060), some subsistence of J. regia was observed along with overall decrease. In SSP2-4.5 (2081–2100) scenario, huge increase in suitable environment was predicted as compared to SSP2-4.5 (2041–2060) (Figs. 3B, C and 4B, C). For future SSP5-8.5 (2041–2060) scenario, survival was observed along with eastward migration. Further, in SSP5-8.5 (2081–2100) future scenario, subsistence was predicted with certain decrease in suitable environment area with eastward migration (Figs. 3D, E and 4D, E).

The model for P. wallichiana predicted most of its suitable area in north of KP in current climate (Fig. 3A). The binary map (Fig. 4A) showed suitable areas in north and east of KP. In SSP2-4.5 (2041–2060) scenario, suitable area increased as compared to current climatic conditions. In SSP2-4.5 (2081–2100) scenario, a similar trend was observed as in current climatic conditions however, binary map (Fig. 4C) predicted the survival of P. wallichiana in existing habitat. For future SSP5-8.5 (2041–2060) scenario, overall increase was observed as compared to current climatic conditions. In SSP5-8.5 (2081–2100) scenario, slight decrease in overall suitable area was observed as compared to SSP5-8.5 (2041–2060) scenario (Figs. 3D, E and 4D, E).

The model for E. camaldulensis predicted most of its suitable area in center of KP (Fig. 3A). According to binary projection (Fig. 4A), other suitable areas were observed in central east of KP. For future SSP2-4.5 (2041–2060) scenario similar areas were projected as in current climatic condition show ever, slight increase in area can be observed. According to binary projection, there is an overall expansion in suitable area for E. camaldulensis is along with migration towards north-east. For SSP2-4.5 (2081–2100) scenario, similar trend was predicted as in SSP2-4.5 (2041–2060) scenario, however, slight expansion in existing suitable areas towards north-east was noticeable (Figs. 3B, C and 4B, C). In SSP5-8.5 (2041–2060) future scenario, suitable areas were considerably changed in the western side of KP. Loss of suitable areas was observed for E. camaldulensis with migration towards east. In SSP5-8.5 (2081–2100) scenario, suitable area in south of KP appeared.

The model for S. modesta under current climatic conditions predicted its potential distribution in central partsof KP (Fig. 3A). Binary projection (Fig. 4A) showed environmental suitability of S. modesta in central and southeastern side of KP. For future scenario SSP2-4.5 (2041–2060), projected results show overall expansion in environmentally suitable area of S. modesta along with migration towards north. In SSP2-4.5 (2081–2100) scenario, a drastic reduction is projected in central and eastern areas as compared to SSP2-4.5 (2041–2060) (Figs. 3B, C and 4B, C). In future SSP5-8.5 (2041–2060) scenario, survival in existing habitat was projected along with migration towards north. In SSP5-8.5 (2081–2100) similar trend was projected as in SSP5-8.5 (2041–2060) (Figs. 3D, E and 4D, E).

The model for P. ciliata (Fig. 3A), under current climate conditions predicted north of KP as the most suitable area. In binary projection (Fig. 4A), suitable areas were predicted in all northern parts of KP. In SSP2-4.5 (2041–2060) scenario (Fig. 3B), suitable areas decreased for P. ciliata in northern parts, however, overall expansion in environmentally suitable areas was observed. In SSP2-4.5 (2081–2100) scenario, suitable habitat reduction for P. ciliata in west along with migration towards east was observed. In SSP5-8.5 (2041–2060) (Fig. 3D), expansions in south were observed along with increase in north. In future SSP5-8.5 (2081–2100) scenario (Fig. 3E), results were similar as in previous timescale SSP5-8.5 (2041–2060) (Fig. 3D). In binary projection (Fig. 4E), some reduction in environmental suitable area was observed.

The model for V. nilotica, under current projection predicted most of its suitable area in south of KP. The binary projection (Fig. 4A) predicted suitable areas from center of KP to south. In SSP2-4.5 (2041–2060) scenario, survival of V. nilotica in existing habitat was observed along with reduction in central area of KP. In SSP2-4.5 (2081–2100) scenario, suitability areas reduced in central and western side of KP. In SSP5-8.5 (2041–2060) future scenario, regional migration has been observed along with reduction in suitable area in south. In SSP5-8.5 (2081–2100) scenario (Fig. 3E), similar trend has been projected as in previous timescale however, environmental suitable areas are increased as compared to previous timescale of SSP5-8.5 (2041–2060) (Figs. 3D, E and 4D, E).

Species at risk of losing habitat

The overall suitable area decreased for C. deodara (Fig. 5) however, C. deodara showed most suitable area under current climate. Dalbergia sissoo showed increase in both future scenarios. Juglans regia exhibited rapid decrease in SSP2-4.5 (2041–2060) scenario, at the end of century under SSP2-4.5 (2081–2100) scenario; it showed the maximum environmentally suitable area. Pinus wallichiana showed overall increase in future scenarios. Eucalyptus camaldulensis showed overall increase in its environmentally suitable area with most suitable area under SSP5-8.5 (2081–2100) scenario and its least suitable area SSP5-8.5 (2041–2060). Senegalia modesta showed in both time periods of SSP2-4.5 and SSP5-8.5 scenarios. Populus ciliate showed its least environmentally suitable area under SSP2-4.5 (2081–2100) scenario however, suitable area increased in the other scenario. Vachellia nilotica showed overall decreasing trend and least suitable area was predicted under SSP5-8.5 (2041–2060).

Figure 5. Suitable area percentage for different species under different scenarios. Few of the scenarios for some species like Juglans regia, Eucalyptus camaludensis, Populas ciliate and Vachellia nilotica show significant increase in suitable area.

Afforestation success

Afforested sites of C. deodara (Fig. 4A), can be seen covering only small part of the suitable area. Few afforestation sites were predicted suitable under SSP2-4.5 (2041–2060) (Fig. 4B). As a result of the reduction of suitable area, a few afforestation sites which were previously predicted as suitable area (Fig. 4B) were simulated to be in non-suitable area under SSP2-4.5 (2081–2100). Most of the afforested sites of D. sissoo fell in non-suitable areas under current climatic conditions and the future scenarios. Juglans regia’s afforested site were predicted non suitable under current climatic conditions and all future scenarios. All afforestation sites of P. wallichiana under current climatic conditions (Fig. 4A) are within suitable area, however all afforested sites predicted to be in the range of non- suitable area under SSP2-4.5 (2041–2060) (Fig. 4B). Further in SSP2-4.5 (2081–2100) scenario few afforested sites of P. wallichiana again appeared in predicted environmentally suitable area. Almost every afforested site was predicted as non-suitable area under both timescales of SSP5-8.5. For E. camaldulensis, few of the afforested sites in central KP were predicted as suitable sites under current climate and future scenarios. Afforested sites of S. modesta under current climatic conditions and SSP2-4.5 (2041–2060) (Fig. 4A, B), were somewhat within predicted suitable areas however, in other three future scenarios few sites were predicted outside the suitable areas. Afforested sites of P. ciliata did not match with the predicted suitable areas under current climate conditions and the two future scenarios. Afforested sites of V. nilotica in the central regions of KP fell in non-suitable area under current climate conditions and both of the future scenarios however, southern afforested sites were predicted within environmentally suitable areas.

Cedrus deodara

Cedrus deodara is considered as landscaping tree and introduced in many parts of world however, the native range of this species in the high mountains of Himalaya where it has established naturally within elevation range of 3,500 to 12,000 feet. The increased demand for the commercial use of this species has threatened in the survival of this species in many parts of its native range (Siddiqui et al. 2013). This species however have been successfully introduced outside of its native range under similar climatic conditions including United State of America. Several varieties are available in trade today for C. deodara. This species generally propagates through seeds however, it can also be grafted as well (Bisht and Bhatt 2016). Considering these characteristics and successful establishing outside of its native range as well as taking the assistance from the model results, the potential habitat suggested by model results, better results can be achieved. Cedrus deodara is also national tree of Pakistan showed decreased in suitable area under scenario SSP2-4.5 (2041–2060) as compared to current climatic conditions and it increased under SSP2-4.5 (2081–2100). Similarly, in SSP5-8.5 (2041–2060) suitable area for C. deodara decreased from current climatic conditions and it increased under SSP5-8.5 (2081–2100) scenario as compared to SSP5-8.5 (2041–2060). Overall decrease in suitable area in all projections were observed, majority of the afforested sites were in non-suitable areas because of slight shift of suitable habitat towards north. Sheikh et al. (2021) also predicted a slight shift of C. deodara towards higher regions due to increase in temperatures.

Dalbergia sissoo

The preferred habitat for this species is the plain area along roadsides, canals, field and forest. This species being the commercially important tree species is facing threats due to conversion of forest, uncontrolled cutting and competition by alien invasive species (Mahatara et al. 2021).

Having the ability to thrive in a variety of soil types from dry to wet but not tolerant to salt. This species has the ability to grow well in riverine environment and mostly prefers alluvial soil, natural sand or gravel conditions. This species have largely been used for restoration of degraded soil because of its ability to improve the soil quality due to its physical, chemical, biological and enzymatic properties. Moreover, soil moisture and water holding capacity for this species is significantly higher (Shah et al. 2010). Considering these characteristics and coupling this with potential habitat predicted by our model result, can help afforestation plan in future.

Dalbergia sissoo, known as “Shisham” tree, is planted for its quality timber in tropical regions (Sangha and Jalota 2005). However, due to bad management practices, inadequate planning, as well as diseases are the major causes of Shisham disappearance from subcontinent (Mukhtar et al. 2015). Suitable area, for D. sissoo under scenario SSP2-4.5 (2041–2060) increased as compared to current climatic conditions, however under scenario SSP2-4.5 (2081–2100) it decreased. Similarly, in SSP5-8.5 (2041–2060) suitable area for D. sissoo increased from current climatic conditions and slightly decreased under SSP5-8.5 (2081–2100) from SSP5-8.5 (2041–2060). Stability of afforested sites of D. sissoo remained critical in all projections.

Juglans regia

Juglans regia is considered an important species for environmental management because of its useful wood, highly nutritious fruit. The species having wide tolerance range for precipitation, temperature and pH is also resistant to cold and temperature down to about –27°C (Sharma et al. 2022). This species requires deep well drained loam and sunny position sheltered from strong wind but also prefers slightly alkaline heavy loam but succeeds in most soils. This species has invasive nature as well which facilitated its range shift and expansions (Sharma et al. 2022). Studies in Europe predicted the northward shift for this species while the distribution is lost in the southern edge under climate change scenarios. This species has native rang in the Himalaya however it has been introduced in many regions of the world indicating successful afforestation of this species. Juglans regia known as Persian walnut, is a valuable woody oil plant, cultivated worldwide for its wood and nuts (Gupta et al. 2019). Our results show that the suitable areas increased according to both future SSPs as compared to current scenario. A decrease in suitable area of J. regia in Europe was predicted by Paź-Dyderska et al. (2021) for three different GCMs, under 3 different climatic projection for (2061–2080) and reported a shift towards northern areas of Europe.

According to the data provided by afforestation project in KP, BTAP, J. regia is afforested on single site and it was not planted in highly suitable area according to our projections.

Pinus wallichiana

Pinus wallichina is native to Bhutan, Nepal, India, and Pakistan. In Pakistan it has been divided into 2 varieties: Var. wallichiana isolated in the moist temperate zone of Murree-Galiat and Azad Kashmir. Artificial plantations are being raised in Kaghan and swat. Var. karakorama isolated in the dry temperate zone of the Northern Areas, Takhte-Sulaiman, Swat, Dir, Chitral, Tirah, and Kurram Plantations are being raised in Kaghan and Swat (Shah et al. 2009). A moderately intolerant tree that grows on a variety of soils but does best on fertile well drained sandy clay to sandy clay loams. It is adapted to a precipitation zone of 300 to 1,500 mm/yr in a temperature range of –20 to 35°C. It prefers a humid cool temperate/arid cold temperate climate at elevations between 1,200 and 3,700 m (Rahman et al. 2020; Shah et al. 2009).

Pinus wallichiana is largely a cold-hardy plant and can establish around –10°C. It prefers well-drained sandy or gravelly loam (Rahman et al. 2020). The species establishes well in position sheltered from strong winds however this species grows fast but it is short-lived in cultivation. Suitable area for P. wallichiana under scenario SSP2-4.5 (2041–2060) increased from current climatic conditions however, it further increased under in rest of the scenarios. Another maxent based study in Himalaya and Hindu-Kush designed on A2a scenario along with GCM. Hadley Centre coupled model version 3, predicted an extensive suitability in current climatic conditions and huge decline in future scenario along with shift towards eastern parts (Ali 2015).

Eucalyptus camaldulensis

Eucalyptus camaldulensiscan be seen as an iconic tree of superlatives. It is the eucalypt with the widest native range, and one of the most widely planted eucalypts around the globe. In South Africa, it is the most widespread and the most aggressively invasive eucalypt (Hirsch et al. 2020). It has many uses, but also causes major impacts. However, little is known about key aspects of its ecology in South Africa, including its invasion history, invasion processes and dynamics, and people’s perceptions of its positive and negative effects on ecosystems. Such knowledge is crucial for developing robust and defendable guidelines for sustainable management of the species. Its adaptation to a wide range of environments has contributed to it becoming one of the most widely cultivated eucalypts across a range of arid, temperate and tropical countries (McDonald et al. 2009). Eucalyptus camaldulensis have been vastly planted in BTAP in saline and water-logged conditions in KP mostly in southern areas (World Wide Fund for Nature Pakistan 2017) though our results suggest suitable areas in central KP. Suitable areas for E. camaldulensis increased in both SSP2s 4.5 and 8.5 which indicates expansion of the species in future with slight shift towards north.

Senegalia modesta

It is a crucial environmental weed in the south-eastern parts of Australia and South Africa. It is recorded as invasive according to the Global Invasive Species Database. It is a very adaptable species that can tolerate relatively dry, low nutrient soils and thrive in better conditions (Hassan and Hamdy 2021). It is planted widely for a range of purposes (e.g., windbreaks, sand dune fixation, fuel woods, fodder production) in West Asia, South America, North and South Africa and parts of the Mediterranean. In Egypt, it was cultivated as an ornamental tree in historical and many public gardens and streets. Recently, 1 million seedlings have been planted along the Mediterranean coast of Egypt for a range of rehabilitation (Hassan and Hamdy 2021; Świerszcz et al. 2023). It escaped from cultivation, spread invasively, and strongly affected biodiversity. Currently, it is utilised in forest plantations in Gebel Elba natural forests, for sand dune fixation, as a bio-fertiliser and windbreaks. Senegalia modesta suitability has shown increase for both SSPs as compared to current scenario. Ali and Begum (2015) found in their study that S. modesta, which is also known as A. modesta, will have large suitable areas in 2080 particularly in Swat valley. Their results are similar to our (2041–2060) timescale predictions but conflicts with (2081–2100) timescale results. Our results predicted that majority of P. ciliata afforested sites were falling under non suitable areas.

Vachellia nilotica

Vachellia nilotica has been planted widely for generations in arid and semi-arid regions of India and Pakistan. It is characteristic of dry regions and does not grow in areas receiving rainfall in excess of about 1,250 mm, or in localities susceptible to frost and cold (Puri et al. 1994). This species is tolerant of salts, arid environments (Minhas et al. 1997), grazing, drought and fire (Healey et al. 1997). Acacias often form the main part of the woody vegetation in the semiarid areas on African savannas (Tybirk 1989). According to our results V. nilotica have maximum suitable area under current climatic conditions and it will decrease in futuristic climate in KP. According to our results, many afforested sites of V. nilotica were predicted under non suitable areas. Taylor et al. (2018) using maxent, predicted that soil, elevation and annual mean temperature as important variables for V. nilotica and future temperature rises may have negative consequences on the distribution of V. nilotica, particularly when winter min-temperatures rise.

Populus ciliata

Is propagated clonally by root suckers but natural regerneration uasually takes place by seeds. This aspect alone cannot be attributed to the prepondaerance of males or females. It was oberved that the male trees were growing in competition with P. wallichiana, C. deodara, P. smithiana trees in the forests and with alnus species. While the females trees are grown in open areas (Siddiqui et al. 1986). This indicates that male trees face more competitions for nutrients as compared with females. It has been observed that thick forest of P. wallichiana, P. smithiana, and C. deodara have little numbers of populus as intermixed vegetation. These types of models can be useful in demonstrating how species can react to restoration efforts (Beatty et al. 2018). Identifying the suitability of the landscape for plantation has always been one of the main challenges for non-forest restoration sites. Identifying areas within the landscape where multiple functions based on various dimensions can be boosted through identifying the appropriateness of the terrain for plantation (Schulz and Schröder 2017). The outputs from most of the studies that used habitat suitability modelling tools or SDM were utilized in biodiversity assessments, foreseeing the consequences of climate change on ecology, indicating potential areas for conservation, habitat restorationandspecies migration (Araújo et al. 2005). To minimize budgetary losses and governmental resources, it is critical that robust research and development institutions, as well as educational institutions are engaged throughout the restoration projects. Our findings, in coming years, can help reforestation programs by indicating where these important species might grow in the future, and which afforested sites needs post plantation care for sustenance.

Our Maxent based results conclude that, changing climate will have very multifaceted effects on certain tree species in future. Afforestation projects without taking into account future predictions, might have a good sustenance rate for current conditions but future fate of these species remains uncertain. Future climatic data integration may offer extra substantiation about the vulnerability of tree species in future and guidelines for conservation of species under future climate scenarios. Selection of tree species for afforestation programs can become a major challenge in coming years as a co-occurrence of tree species for resource utilization with changing climate cannot be comprehended without spatial techniques. Maxent based studies are mostly used for predicting the current and future habitats of tree species but the study also indicates that the SDMs can be useful for inspecting the fate of already afforested programs.

SSP: Shared socio-economic pathway

SDM: Species distribution model

BTAP: Billion Tree Tsunami Afforestation Project

KP: Khyber Pakhtunkhwa

GCM: Global climate models

CMIP6: Coupled Model Inter comparison Project

DEM: Digital elevation model

AUC: Area under the curve

MAD carried out the research, RR helped with running the model, MS helped with data collection, MD helped with valuable input and proofreading the manuscript and SS conceived the research idea and supervised the research.

National University of Sciences & Technology (NUST) student research grant.

Not applicable.

Not applicable.

Not applicable.

The authors declare that they have no competing interests.