Predicting the Potential Geographic Distribution and Habitat Suitability of Two Economic Forest Trees on the Loess Plateau, China

The Loess Plateau is one of the most fragile ecosystems in the world. In order to increase the biodiversity in the area, develop sustainable agriculture and increase the income of the local people, we simulated the potential geographic distribution of two economic forest trees (Malus pumila Mill and Prunus armeniaca L.) in the present and future under two climate scenarios, using the maximum entropy model. In this study, the importance and contributions of environmental variables, areas of suitable habitats, changes in habitat suitability, the direction and distance of habitat range shifts, the change ratios for habitat area and land use proportions, were measured. According to our results, bioclimatic variables, topographic variables and soil variables play a significant role in defining the distribution of M. pumila and P. armeniaca. The min temperature of coldest month (bio6) was the most important environmental variable for the distribution of the two economic forest trees. The second most important factors for M. pumila and P. armeniaca were, respectively, the elevation and precipitation of the driest quarter (bio17). At the time of the study, the area of above moderately suitable habitats (AMSH) was 8.7967 × 104 km2 and 11.4631 × 104 km2 for M. pumila and P. armeniaca. The effect of Shared Socioeconomic Pathway (SSP) 5-85 was more dramatic than that of SSP1-26. Between now and the 2090s (SSP 5-85), the AMSH area of M. pumila is expected to decrease to 7.5957 × 104 km2, while that of P. armeniaca will increase to 34.6465 × 104 km2. The suitability of M. pumila decreased dramatically in the south and southeast regions of the Loess Plateau, increased in the middle and west and resulted in a shift in distance in the range of 78.61~190.63 km to the northwest, while P. armeniaca shifted to the northwest by 64.77~139.85 km. This study provides information for future policymaking regarding economic forest trees in the Loess Plateau.


Introduction
The global ecosystem is under great pressure resulting from both climate change and social development [1]. A number of important studies have proven that climate change is having a profound impact on ecosystem biodiversity, at both global and local scales [2]. The global average temperature has risen by 0.1 • C per decade since the 1950s and is projected to increase by 1.5-2.1 • C by 2050, with a maximum rise of 2.6-4.8 • C by the end of this century [3,4]. Climate change (i.e., temperature, precipitation and radiation) leads to land use and land cover change (LUCC), which in turn leads to changes in the climate, particularly in the form of temperature changes. Both climate change and LUCC have significant impacts on the distribution ranges and spatial living patterns of biological organisms.
The importance and contributions of climate variables in explaining the distribution of species have long been recognized [5]. In the past few decades, with technological

Species Occurrence Data
In order to obtain the occurrence records for the northern region of the Yangtze River in China (except Qinghai, Tibet and Xinjiang), we applied the location data from the Chinese Virtual Herbarium (a platform supported by the Chinese government for sharing the digital information of nearly 8 million plant specimens from more than 100 domestic herbarium, (https://www.cvh.ac.cn, accessed on 23 November 2020) and retrieved 501 records (260 of M. pumila and 241 of P. armeniaca). Since the prediction of plant species' geographic distribution requires high-accuracy occurrence data, we first excluded duplicate records [4] and those records whose locations were above the county level. Then, we used the Google Earth map of Earth Online (https://www.earthol.com/, accessed on 29 December 2020) to determine the latitude and longitude coordinates according to the position description of the specimen, particularly when the occurrence records lacked exact coordinates. To reduce spatial autocorrelation and improve the accuracy of model prediction, we used the SDM toolbox Ver 2.4 of ArcGIS 10.2 (ESRI, Redlands, CA, USA) to rarefy occurrence data for SDMs. The rarefication ensures that the distance between any two occurrence records is greater than 5 km. Finally, we obtained 176 geographic locations of M. pumila and P. armeniaca (Figure 1) for the Maxent model.

Species Occurrence Data
In order to obtain the occurrence records for the northern region of the Yangtze Rive in China (except Qinghai, Tibet and Xinjiang), we applied the location data from the Ch nese Virtual Herbarium (a platform supported by the Chinese government for sharing th digital information of nearly 8 million plant specimens from more than 100 domestic her barium, (https://www.cvh.ac.cn, accessed on 23 November 2020) and retrieved 501 rec ords (260 of M. pumila and 241 of P. armeniaca). Since the prediction of plant species' geo graphic distribution requires high-accuracy occurrence data, we first excluded duplicat records [4] and those records whose locations were above the county level. Then, we use the Google Earth map of Earth Online (https://www.earthol.com/, accessed on 29 Decem ber 2020) to determine the latitude and longitude coordinates according to the positio description of the specimen, particularly when the occurrence records lacked exact coor dinates. To reduce spatial autocorrelation and improve the accuracy of model prediction we used the SDM toolbox Ver 2.4 of ArcGIS 10.2 (ESRI, Redlands, CA, USA) to raref occurrence data for SDMs. The rarefication ensures that the distance between any tw occurrence records is greater than 5 km. Finally, we obtained 176 geographic locations o M. pumila and P. armeniaca (Figure 1) for the Maxent model.

Environmental Variables
To describe future climate change [30], we used the same bioclimatic variable projec tion (2.5 arc-min resolution) of the Beijing Climate Center Climate System Model (BCC CSM2-MR) in phase 6 of the Coupled Model Intercomparison Projects (CMIP6), whic has been widely used in East Asia [31,32]. Among four emission scenarios of the Share Socioeconomic Pathways (SSP), SSP1-26 intends to limit warming below 1.5 °C, whil SSP5-85 predicts an increase of 5 °C by 2100 [33,34]. Two of the climate scenarios, SSP1-2 (the optimistic case, 2.6 Wm -2 level) and SSP5-85 (the worst case, 8.5 Wm -2 level), showin the two most extreme scenarios of future climate change [4], were finally, selected as th future variable types for the Maxent model [35].

Environmental Variables
To describe future climate change [30], we used the same bioclimatic variable projection (2.5 arc-min resolution) of the Beijing Climate Center Climate System Model (BCC-CSM2-MR) in phase 6 of the Coupled Model Intercomparison Projects (CMIP6), which has been widely used in East Asia [31,32]. Among four emission scenarios of the Shared Socioeconomic Pathways (SSP), SSP1-26 intends to limit warming below 1.5 • C, while SSP5-85 predicts an increase of 5 • C by 2100 [33,34]. Two of the climate scenarios, SSP1-26 (the optimistic case, 2.6 Wm −2 level) and SSP5-85 (the worst case, 8.5 Wm −2 level), showing the two most extreme scenarios of future climate change [4], were finally, selected as the future variable types for the Maxent model [35].
We initially selected 24 environmental factors that may affect the geographic distribution of M. pumila and P. armeniaca to model the present and future species distribution patterns. These factors included 19 bioclimatic variables (bio1-bio19) with 2.5 arc-min spatial resolution (~5 km) from WorldClim v2.1 (https://www.worldclim.org, accessed on 18 May 2020) [36], one digital elevation model datum (DEM,~1 km) and three soil texture data (~1 km) from RESDC (http://www.resdc.cn/, accessed on 19 May 2020) and one soil type datum from FAO (http://www.fao.org/soils--portal/en, accessed on 11 May 2020). With the help of the toolbox in ArcGIS, three extra topographic factors (slope, aspect, curvature) were obtained by using DEM data (elevation). To avoid over-fitting [29,37], 8 out of 27 variables were selected (more details are given in [38]). Based on the physiology requirements [30] of M. pumila and P. armeniaca, we selected four additional bioclimatic variables. Finally, 12 vital environmental variables were obtained for the Maxent model, including six bioclimatic variables, four topographic variables and two soil variables (Table 1). To describe the land use situation, we chose the observational data (300 m spatial resolution) of 2010 from Copernicus (https://climate.copernicus.eu/, accessed on 7 August 2020). As shown in Table S1, we combined and reclassified land use types with similar functions and finally obtained seven types: cropland, grassland, forest, mosaic area, bare, urban and built-up and water.

Model Implementation and Evaluation
In this study, we used Maxent 3.4.1 [6] to predict the potential geographic distribution of species. To ensure the accuracy of Maxent, with the help of ArcGIS, we set the environmental variables to the same resolution (2.5 arc-min) and the same projected coordinate system (Asia North Alberts Equal Area Conic, WGS1984). In addition to the recommended default parameters [6], we ran the model with 30% test data (70% of training data) and 5 replicates. We firstly use maximizing of the sum of sensitivity and specificity (Maxsss) to delineate the binary presence/absence maps [39]. The presence distribution maps were then reclassified and divided into lowly suitable, moderately suitable and highly suitable areas by breakpoints 0.4 and 0.6.

Spatial and Statistical Analysis
With the help of the raster calculator, we calculated the changes in suitability. In the process of describing the increase and decrease maps of M. pumila and P. armeniaca, we define the interval period as "Pt-Pp", while Pt represents the target period and Pp represents the previous period. The reclassified suitability change was divided into seven conditions: high decrease, middle decrease, low decrease, no change, low increase, middle increase and high increase. Specifically, the no change parts were divided into two specific In order to determine the moderately suitable and highly suitable area proportions of land use in different periods, we overlaid the reclassified suitable maps and the LUCC maps and finally obtained the suitability conditions with land use in three periods. Since it is difficult to obtain relatively accurate simulation results of LUCC over a long time span, we assumed that the land use will not change in the next 100 years. By considering the gravity of suitability, we used the mean center tool to calculate the distribution maps of M. pumila and P. armeniaca under the climate scenarios and finally obtained the centroids in different periods. The change ratio (Cr) represents the suitable habitat area change from the present (Ap) to the target period (At) and it can be calculated using the following equation: Cr = (At − Ap)/Ap [30]. The combination of periods (2050s and 2090s) and climate scenarios (SSP1-26 and SSP5-85) represents the condition of the period under the corresponding climate scenario.

The Contribution and Importance of Environmental Variables
The accuracy of Maxent was evaluated based on the AUC values. During the Maxent simulation, all the AUC values were higher than 0.9. The average AUC values were 0.946 ± 0.02 and 0.959 ± 0.02 (mean ± SD) for M. pumila and P. armeniaca, indicating that Maxent displayed excellent prediction accuracy in our simulations. These results also indicated that the bioclimatic data, topographic data and soil data played an important role in developing the model. Among all environmental variables, bio6 and elevation played important roles in determining the future distribution of M. pumila in both the 2050s and 2090s. The contribution of bio6 and elevation exceeded 72%, while their importance exceeded 73% in all simulations of the 2050s and 2090s under SSP1-26 and SSP5-85 ( Figure 2A,B). However, it should be noted that both bio6 and elevation had a wide range of contribution rate and importance rate values in the SSP1-26 scenario (i.e., the contribution of bio6 ranged from 33.8% to 48.4% in the 2050s, the importance of bio6 varied from 29.1% to 51.4% in the 2050s, the importance of elevation ranged from 22.2% to 49.2% in the 2050s and the importance of elevation ranged from 19.0% to 43.0% in the 2090s). Interestingly, compared with SSP1-26, the contribution rate and importance rate of bio6 and elevation in the SSP5-85 scenario were relatively stable, particularly the contribution rate of bio6 (43.4-46.7%) ( Figure 2B).
During the simulation of P. armeniaca, the three highest-ranked factors were bio6, bio 12 and bio17, whose cumulative relative contribution and importance exceeded 65% and 76%, respectively. The average contribution rate of bio6 was the highest in both the 2050s and 2090s under the SSP1-26 and SSP5-85 scenarios. Among all environmental variables used to predict the distribution of P. armeniaca, bio12 ranked third, but it had the largest range (i.e., the importance rate varied from 2.2% to 30.3% under the SSP5-85 climate scenario) in both the 2050s and 2090s.
For M. pumila and P. armeniaca, bioclimatic variables are the most important variables affecting their growth and development, followed by topographic variables and soil variables. Among all environmental factors, the contribution rate and importance rate of bio1, bio5 and aspect were close to 0. For M. pumila, the contribution rate of topographic and soil factors showed the following order: elevation > sand > slope >curvature > soil > aspect. However, the contribution ranking of topographic and soil factors for P. armeniaca was as follows: elevation > slope > curvature > soil > sand > aspect.

Suitable Habitats for M. Pumila and P. Armeniaca
We obtained the distribution maps of M. pumila and P. armeniaca in the present and future (SSP1-26 and SSP5-85 scenarios, Figure 3). The suitability was reclassified into four categories (unsuitable, low suitable, moderately suitable and highly suitable) and the corresponding habitats were classified as unsuitable habitat (UH), low suitable habitat (LSH), moderately suitable habitat (MSH) and highly suitable habitat (HSH).

Suitable Habitats for M. pumila and P. armeniaca
We obtained the distribution maps of M. pumila and P. armeniaca in the present and future (SSP1-26 and SSP5-85 scenarios, Figure 3). The suitability was reclassified into four categories (unsuitable, low suitable, moderately suitable and highly suitable) and the corresponding habitats were classified as unsuitable habitat (UH), low suitable habitat (LSH), moderately suitable habitat (MSH) and highly suitable habitat (HSH).
For M. pumila, the area of UH, LSH, MSH and HSH was 43.7759 × 10 4 km 2 , 12.3149 × 10 4 km 2 , 7.9449 × 10 4 km 2 and 0.8518 × 10 4 km 2 , respectively ( Figure 4A, Table S2). The suitable area of M. pumila showed a difference change under SSP1-26 and SSP5-85 on the Loess Plateau. Under SSP1-26, the area of suitable habitats (SH: LSH, MSH and HSH) showed a slight trend of continuous increase from the present period to the end of this century (Table S2); however, the area of above moderately suitable habitats (AMSH: MSH and HSH) dropped from 8.7967 × 10 4 km 2 (present) to 7.8458 × 10 4 km 2 for the 2050s and then increased to 18.3296 × 10 4 km 2 for the 2090s. Interestingly, compared with SSP1-26, the area of SH under SSP5-85 did not maintain a stable and sustained change and the area of SH increased for the 2050s and decreased dramatically for the 2090s. The area of AMSH was consistent with the change trend of SH, which increased from the present to 11.1493 × 10 4 km 2 for the 2050s and then dropped to 7.5957 × 10 4 km 2 for the 2090s.   (Table S2); however, the area of above moderately suitable habitats (AMSH: MSH and HSH) dropped from 8.7967 × 10 4 km 2 (present) to 7.8458 × 10 4 km 2 for the 2050s and then increased to 18.3296 × 10 4 km 2 for the 2090s. Interestingly, compared with SSP1-26, the area of SH under SSP5-85 did not maintain a stable and sustained change and the area of SH increased for the 2050s and decreased dramatically for the 2090s. The area of AMSH was consistent with the change trend of SH, which increased from the present to 11.1493 × 10 4 km 2 for the 2050s and then dropped to 7.5957 × 10 4 km 2 for the 2090s.

Changes in Habitat Suitability and Range Shifts
The change in suitability was divided into seven levels ( Figure 5) based on the reclassified suitability of habitats. For M. pumila, by the 2050s (SSP1-26), the area of habitat suitability was predicted to mainly decrease in the south of the Loess Plateau (a low decrease in the south and a dramatic decrease in the southeast); the suitability remained largely unchanged in the central, northern and southwestern areas (the central area was suitable habitat, while the north and southwest were unsuitable areas) and there was a scattered increase in suitability in the central area ( Figure 5). By the 2090s (SSP1-26), the area of habitat suitability was predicted to decrease in the south and southeast of the Loess Plateau and increase in the central and western areas. Specifically, compared with SSP1-26, the habitat suitability change in M. pumila showed different trends under SSP5-85. By the 2050s, in the south of the Loess Plateau, a dramatic decrease in habitat suitability was predicted over a larger area, while the central and western areas showed a dramatic increase ( Figure 5). However, by the 2090s, the area of suitability was predicted to dramatically decrease in the south of the Loess Plateau, while the area of suitability will dramati-   Table S2), while under SSP5-85, the area in the two periods was 26.6638 × 10 4 km 2 and 34.6465 × 10 4 km 2 . The area of MSH and HSH increased greatly; in particular, the area of HSH under the SSP5-85 scenario increased by 550% compared to the present area (Table S2).

Changes in Habitat Suitability and Range Shifts
The change in suitability was divided into seven levels ( Figure 5) based on the reclassified suitability of habitats. For M. pumila, by the 2050s (SSP1-26), the area of habitat suitability was predicted to mainly decrease in the south of the Loess Plateau (a low decrease in the south and a dramatic decrease in the southeast); the suitability remained largely unchanged in the central, northern and southwestern areas (the central area was suitable habitat, while the north and southwest were unsuitable areas) and there was a scattered increase in suitability in the central area ( Figure 5). By the 2090s (SSP1-26), the area of habitat suitability was predicted to decrease in the south and southeast of the Loess Plateau and increase in the central and western areas. Specifically, compared with SSP1-26, the habitat suitability change in M. pumila showed different trends under SSP5-85. By the 2050s, in the south of the Loess Plateau, a dramatic decrease in habitat suitability was predicted over a larger area, while the central and western areas showed a dramatic increase ( Figure 5). However, by the 2090s, the area of suitability was predicted to dramatically decrease in the south of the Loess Plateau, while the area of suitability will dramatically increase in the southwest of the Loess Plateau. From the 2050s to the 2090s, in the central and southern areas of the Loess Plateau, M. pumila was predicted to show a slight increase in habitat suitability under SSP1-26 and a slight decrease under SSP5-85. For the habitat suitability change in P. armeniaca, from the southeast to northwest, the intensity of the change in suitability varied from a high decrease to a high increase in all cases. The suitability of P. armeniaca remained constant in the middle of the Loess Plateau and dramatically increased in the northern and northwestern regions. SSP5-85 displayed a change of greater intensity than SSP1-26.

the Dominant Composition and Proportion of Land Cover
We obtained the preferable land cover proportions of M. pumila and P. armeniaca under SSP1-26 and SSP5-85 (Figure 7). At present, for M. pumila, the dominant land cover types of AMSH were, in order, cropland, forest and mosaic areas (Figure 7). Under SSP1-26, the proportion of cropland decreased from 34.68% for the 2050s to 32.51% for the 2090s and the proportion of grassland increased from 23.28% for the 2050s to 28.42% for the 2090s. SSP5-85 showed a stronger change range than SSP1-26. Compared with the proportion of the present, forest type showed an initial increasing trend and then decreased under both the SSP1-26 and the SSP5-85 scenarios. For P. armeniaca, the dominant land cover types were, in order, cropland, mosaic area and forest at present (Figure 7). The proportions of grassland showed an increasing trend from the present to the 2050s and the 2090s and SSP5-85 saw a greater increase in proportion than SSP1-26. By the 2090s, the proportion of cropland was predicted to decrease to 30.35% under SSP1-26 from 39.52% at present and to decrease to 21.78% under SSP5-85; the proportion of mosaic area showed a stable trend, which was predicted to change from 22.38% at present to 22.95% (2050s) and 24.08% (2090s) under SSP1-26 and to 22.89% (2050s) and 22.26% (2090s) under SSP5-85 (Figure 7). and SSP5-85 saw a greater increase in proportion than SSP1-26. By the 2090s, the proportion of cropland was predicted to decrease to 30.35% under SSP1-26 from 39.52% at present and to decrease to 21.78% under SSP5-85; the proportion of mosaic area showed a stable trend, which was predicted to change from 22.38% at present to 22.95% (2050s) and 24.08% (2090s) under SSP1-26 and to 22.89% (2050s) and 22.26% (2090s) under SSP5-85 ( Figure 7).

Environmental Factors Shaping the Distribution of Species
We simulated the geographic distribution of M. pumila and P. armeniaca at present and in the 2050s and 2090s by Maxent under two emission scenarios (SSP1-26, SSP5-85). Our results, which indicate robust and accurate model performance (AUC > 0.9), reveal that bioclimatic factors, together with topographic factors and soil factors, play an important role in affecting species distribution, which is consistent with previous studies [7,40]. Among the various environmental variables, min temperature of coldest month (bio6) and elevation are considered the most obvious limitations for the distribution of M. pumila, while for P. armeniaca, the limiting factors are the min temperature of coldest month (bio6) and precipitation of the driest quarter (bio17). Previous studies have proven

Environmental Factors Shaping the Distribution of Species
We simulated the geographic distribution of M. pumila and P. armeniaca at present and in the 2050s and 2090s by Maxent under two emission scenarios (SSP1-26, SSP5-85). Our results, which indicate robust and accurate model performance (AUC > 0.9), reveal that bioclimatic factors, together with topographic factors and soil factors, play an important role in affecting species distribution, which is consistent with previous studies [7,40]. Among the various environmental variables, min temperature of coldest month (bio6) and elevation are considered the most obvious limitations for the distribution of M. pumila, while for P. armeniaca, the limiting factors are the min temperature of coldest month (bio6) and precipitation of the driest quarter (bio17). Previous studies have proven that winter temperature, especially the minimum temperature, has huge effect on the distribution of species, especially for plants at high altitudes [41]. Temperature increase in the future may change the distribution patterns of M. pumila and P. armeniaca [42]. As economic trees that have been widely cultivated in the Loess Plateau, although M. pumila and P. armeniaca are limited by similar bioclimatic factors, they are restricted by different topographic and soil factors. Our simulation shows that, by the 2050s, compared with SSP5-85 (9.33%), the area where the suitability of M. pumila remained unchanged was predicted to maintain a relatively higher percentage (17.67%, Figure 5) under the SSP1-26 scenario, which indicates that M. pumila are sensitive to climate change. This is consistant with the previous studies in Northern India, which have shown that the elevation of apple orchards risen from 1200-1500 m to 1500-2500 m in two decades [41]. With the high contribution and importance of elevation and the low contribution and importance of bio12, bio16 and bio17, M. pumila exhibits much greater tolerance to precipitation and higher sensitivity to the cultivation altitude [42,43]. The high contribution and importance of bio17 indicates that P. armeniaca is more sensitive to precipitation. Therefore, P. armeniaca can be used to occupy the original position of M. pumila [27,42] and future cultivation may require a variety of irrigation to compensate [41]. The texture and type of soil is essential to support plants and provide nutrients. The difference in the contributions and importance of sand (soil texture) and soil (soil type) reveals that, for M. pumila, sand has a greater effect on the geographic distribution than soil, while soil plays a relatively significant role in affecting the distribution of P. armeniaca.

The Distribution of M. pumila and P. armeniaca
The area of AMSH showed differences in the change trend between M. pumila and P. armeniaca under the two scenarios (SSP1-26, SSP5-85). The AMSH area of M. pumila showed the opposite trend under the two scenarios, while that of P. armeniaca showed an increasing trend over time (Figure 4, Table S2). We speculate that the preference for colder temperatures and the limited tolerance of increases in temperature in M. pumila [41] are possible reasons for the contrasting change trends under the two climate scenarios. For P. armeniaca, the increasing intensity of the predicted AMSH area under SSP5-85 was greater than that under SSP1-26, indicating that P. armeniaca had a higher temperature tolerance than M. pumila. Considering the future climate conditions, P. armeniaca would be more suitable for planting in the vast area of the Loess Plateau.
In the biological world, the distribution of species is often affected by many factors simultaneously [1]. It is difficult to quantitatively measure the contributions of various factors to the range of changes in species using general methods and we believe that qualitative analysis is often more important in order to elucidate specific problems. Generally, a smaller temperature increase and higher carbon dioxide concentrations are more conductive to plant growth, but by the 2050s, the temperature accumulation of SSP5-85 may have exceeded the optimal growth temperature for M. pumila, resulting in a decrease in the AMSH area. Recent studies have suggested that temperatures above the optimal temperature range force stomata closure and reduce the overall metabolic activity of plants [44]. In addition, global warming may lead to a redistribution of precipitation across regions, which in turn will significantly alter plant distribution patterns [33].
Under both the SSP1-26 and SSP5-85 scenarios, compared with the centroid of the present, both M. pumila and P. armeniaca shift to the northwest, which is consistent with the direction of the increasing altitude of the Loess Plateau. Under the same climate scenario, M. pumila and P. armeniaca showed the same trend: the range shift distance for the 2090s was larger than that for the 2050s. Many species have been reported to shift poleward in latitude and upward in altitude in many cases, at a large scale, under the continued pressure of global warming [45][46][47]. For climate-sensitive species, similar climate change pressure may result in dramatic expansion or contraction [48], which may result in a longdistance range shift. However, at local or regional scale, species are greatly affected by local factors, which often leads to diverse shifts [45]. Thus, unusual shifts are likely to become explicable when local conditions are taken into account.

Establishment and Conservation Strategies of Economic Forest Trees in the Loess Plateau
It is widely accepted that changes in the global climate will inevitably affect land use, land cover and the distribution of natural species [15,28]. Previous studies have proven that increasing the biodiversity of ecosystems can increase their resistance to climate changes [49]. In the extensive, ecologically sensitive area of the Loess Plateau, policymakers are seeking novel ways to reach the goals of balanced local food security, reduction of poverty and development of sustainable ecosystems [50]. Combining the vegetation restoration and national projects such as the GGP to increase the cultivated area of local economic tree species is a feasible approach to achieve these goals. However, as different species have different levels of tolerance to the environment of the Loess Plateau, it is vital to manage the suitable areas of target species according to their environmental variables in order to increase the possibility of their successful cultivation [18]. In addition, the introduction of some species with good economic benefits can be considered to increase the income of the local people and increase the biodiversity of the Loess Plateau ecosystem. However, when selecting the species to be introduced, the interspecific relationships between the introduced and native species should be taken into consideration.
To ensure that the Loess Plateau's ecosystems can fulfil their regional functions, it is necessary to adopt tailored and comprehensive methods for these fragile ecological areas [51]. Vegetation restoration in different regions may have different purposes. For species with stable environmental supply, such as mangroves in the southeast of China [18], areas of high suitability are the target sites for the establishment of protected areas. In contrast, vegetation restoration in the Loess Plateau aims not only to increase the diversity of the region but also to increase the income of the local people. Until the end of this century, for both M. pumila and P. armeniaca, the land use type with the largest suitable area should be converted from cropland to grassland. However, it should be noted that cropland, grassland, forests and mosaic areas still occupy more than 90% of the AMSH. When considering where to plant economic trees in the Loess Plateau, we should not restrict our selection to the potential areas with moderate and high suitability; in special cases, areas with low suitability could also be appropriate. Combined with the goal of the GGP-to conserve water and soil and increase biodiversity-it is feasible to cultivate economic forest trees in grassland with sufficient moisture, on roadsides in built-up areas and on croplands with steep slopes [24] and patch planting is recommended in forest and mosaic areas.

Model Application and Optimization in Predicting the Distribution of Species
In the present study, the Maxent model was used to predict the distribution of M. pumila and P. armeniaca under SSP1-26 and SSP5-85. Though Maxent showed high accuracy in the statistical methods, there are still some uncertainties. For Maxent, to avoid overfitting, we rarefied the environmental factors. However, those environmental factors that were screened out may contain important information that affects the distribution of species [52]. Moreover, the application of the dynamic model to analyze the influence of individual factors will further improve the accuracy of the simulation. Thus, more rational algorithms are needed. In addition, the extensive collection of occurrence data and high-resolution maps [30] will further improve the robustness of species simulation.
In addition to the common environmental variables [53,54], we tried to incorporate land use and land cover into suitable habitats of economic trees to better understand the impact of climate change and provide a complement for better and more accurate decision-making. However, a long-term LUCC simulation may reduce the accuracy of the simulation, because LUCC is strongly affected by social factors, such as policy, in addition to climate change [28]. It was unclear whether we could simulate LUCC for the next 100 years with relative accuracy, so we used the same LUCC map for all three periods considered in the study. However, in order to obtain high-accuracy LUCC maps, more relevant factors (e.g., the road, elevation and location of the farmland protection zone) must be taken into account to improve the stability of species habitat simulation, which is necessary for future research.

Conclusions
This study investigated the geographical distribution of M. pumila and P. armeniaca on the Loess Plateau using the Maxent model in three periods (present, 2050s and 2090s) under two climate scenarios (SSP1-26 and SSP5-85). Bio6 was considered the dominant factor affecting the distribution of M. pumila and P. armeniaca. M. pumila showed greater sensitivity to elevation and temperature changes. Bio17 played a significant role in the distribution of P. armeniaca. Soil type had a greater effect on the distribution of P. armeniaca, while soil texture had a greater effect on M. pumila. The AMSH area of M. pumila is 8.7967 × 10 4 km 2 at present. By the 2090s, the AMSH area was predicted to increase to 18.3296 × 10 4 km 2 under SSP1-26 and decrease to 7.5957 × 10 4 km 2 under SSP5-85. The SSP5-85 results in a higher range shift in distance than the SSP2-26, for both M. pumila and P. armeniaca. The suitability of M. pumila decreases in the south and southeast and increase in the middle and western areas of the Loess Plateau. Climate change leads to northwest range shifts in M. pumila, with a distance range of 78.61 km to 190.63 km. Compared to the present area