Identifying Potential Planting Sites for Three Non-Native Plants to Be Used for Soil Rehabilitation in the Tula Watershed

: The Tula watershed in Mexico, located in a semiarid and sub-humid climate zone, is experiencing intensive population growth, the expansion of mining concessions for construction materials, and agricultural and urban development, resulting in the degradation of soils and vegetation and a greater demand on natural resources. The aims of this study were to evaluate the survival rates and identify potential habitats within the Tula watershed for planting three non-native forage species ( Atriplex canescens , Cynodon dactylon, and Leucaena collinsii ) using the Kaplan-Meier estimator and the MaxEnt model with the purpose of rehabilitating degraded soils via agroforestry systems. There were 19 edaphoclimatic variables used and the occurrences of three species, obtained from the GBIF, MEXU, and SNIB databases. The models generated with MaxEnt were very accurate (area under the curve [AUC] ≥ 0.7). The species Atriplex canescens and Cynodon dactylon showed areas of potential planting sites (>0.4) and high survival rates (80% and 92%, respectively). The species Leucaena collinsii presented areas with lower potential planting (<0.4) but registered the greater survival rate (100%). The results provide a solid basis to evaluate the survival rates of forage species within potential planting sites in the Tula watershed using agroforestry systems to rehabilitate degraded soils.


Introduction
Soil degradation is the loss of the intrinsic physical, chemical, and/or biological qualities of soil either by natural or anthropogenic processes, which reduces the ecosystem's ability to produce goods and services [1]. The causes of soil degradation include agricultural, industrial, and commercial pollution, the loss of arable land caused by urbanization, overgrazing, and unsustainable agricultural practices, and long-term climatic changes [2]. Soil rehabilitation and conservation can be carried out by vegetative, agronomic, mechanical or a combination of these practices. Vegetative and agronomic practices are commonly employed because of their low costs and numerous benefits. Agroforestry is one of the vegetative practices [3] that can significantly improve degraded soils by reducing surface runoff and erosion, reducing the effect of raindrops on the soil surface, and improving soil quality and recovery [4,5] as well as ensuring the continuity of tree strata [6]. Agroforestry requires that the chosen species (native and non-native) for rehabilitating are not only compatible with the biophysical requirements [7,8] but also with the social, economic, and cultural environment of the proposed planting site [9]. The main goal when selecting species for rehabilitation purposes is to achieve high plant survival and growth rates during the field establishment stage. The adaptive capacity of the species is determined by the The average annual temperature is 14 • C while the temperature of the coldest month varies between −3 • C and 18 • C, and the temperature of the hottest month is less than 22 • C. The predominant vegetation consists mainly of forests of tascate, pine, and oak, and crasicaule scrubs [19,[44][45][46]. There are two main climates within the watershed: semiarid temperate (28% of the total area) and temperate sub-humid (72% of the total area). In the temperate semiarid climate, annual precipitation ranges from 500 to 600 mm per year. The predominant soils are haplic phaeozem and eutrophic cambisol [46,47], both of medium texture (loamy). In the sub-humid temperate climate, annual precipitation ranges from 600 to 1000 mm per year with rainfall in summer and the percentage of winter rainfall is 5 to 10.2% of the annual total. The precipitation in the driest month is less than 40 mm and the predominant soils are medium-textured haplic phaeozem and fine-textured pelic vertisol [19,[44][45][46].

Study Area
The study area was the Tula watershed located in the states of Mexico, Hidalgo, and Tlaxcala ( Figure 1), with elevations from 2333 to 3223 m asl and an area of 1037.66 km 2 . The average annual temperature is 14 °C while the temperature of the coldest month varies between −3 °C and 18 °C, and the temperature of the hottest month is less than 22 °C. The predominant vegetation consists mainly of forests of tascate, pine, and oak, and crasicaule scrubs [19,[44][45][46]. There are two main climates within the watershed: semiarid temperate (28% of the total area) and temperate sub-humid (72% of the total area). In the temperate semiarid climate, annual precipitation ranges from 500 to 600 mm per year. The predominant soils are haplic phaeozem and eutrophic cambisol [46,47], both of medium texture (loamy). In the sub-humid temperate climate, annual precipitation ranges from 600 to 1000 mm per year with rainfall in summer and the percentage of winter rainfall is 5 to 10.2% of the annual total. The precipitation in the driest month is less than 40 mm and the predominant soils are medium-textured haplic phaeozem and fine-textured pelic vertisol [19,[44][45][46].

Experimental Site
The field study of plant survival was established in Ejido San Felipe Teotitlán, municipality of San Felipe Teotitlán, located in the semiarid temperate zone within the Tula watershed, in an area of 1500 m 2 (50 m long × 30 m wide) located at an altitude of 2473 m asl [48]. As the soil in the area is very shallow, holes were dug at 2 m between plants in rows of 10 m and a depth of 30 cm ( Figure 2); before planting, the holes were filled with soil from the same region (Phaeozem) [46,47] to promote root development. During the first week of May 2018, twenty-five plants of L. collinsii (two years old) and seedlings of A. canescens were transplanted into the holes. Stolons of C. dactylon (20 cm long) were transplanted during the first week of July 2019, according to recommendations made by Chaturvedi et al. [49] who suggested planting grasses one year after the establishment of trees and shrubs. The variables plant height, stem diameter, and canopy diameter were measured on the species A. canescens and L. collinsii, during the months of July 2018 and

Experimental Site
The field study of plant survival was established in Ejido San Felipe Teotitlán, municipality of San Felipe Teotitlán, located in the semiarid temperate zone within the Tula watershed, in an area of 1500 m 2 (50 m long × 30 m wide) located at an altitude of 2473 m asl [48]. As the soil in the area is very shallow, holes were dug at 2 m between plants in rows of 10 m and a depth of 30 cm ( Figure 2); before planting, the holes were filled with soil from the same region (Phaeozem) [46,47] to promote root development. During the first week of May 2018, twenty-five plants of L. collinsii (two years old) and seedlings of A. canescens were transplanted into the holes. Stolons of C. dactylon (20 cm long) were transplanted during the first week of July 2019, according to recommendations made by Chaturvedi et al. [49] who suggested planting grasses one year after the establishment of trees and shrubs. The variables plant height, stem diameter, and canopy diameter were measured on the species A. canescens and L. collinsii, during the months of July 2018 and July 2019, as well as their survivorship during the period July 2018 to October 2019. The variables plant height and canopy diameter were measured on the species C. dactylon, in August and October 2019, and its survivorship was evaluated from August to October 2019.
July 2019, as well as their survivorship during the period July 2018 to October 2019. The variables plant height and canopy diameter were measured on the species C. dactylon, in August and October 2019, and its survivorship was evaluated from August to October 2019.

Target Species
The Atriplex canescens (Pursh) Nutt. species known as Chamizo (fourwing saltbush) is a shrub in the Amaranthaceae family native to the arid zones of North America, distributed from Canada to central Mexico [50]. This species is an important food source, both for both domestic animals (mainly cows and goats) and for wild animals [51], and it is used for reforestation, soil conservation, and soil recovery [21]. Cynodon dactylon (L.) Pers. (Bermuda grass) is a perennial grass in the Poaceae family native to tropical Africa, Eurasia, India, and Malaysia and is widely distributed in tropical and subtropical regions in the world [52]. This species is used as a source of food for livestock as well as a phytoremediation tool for salt-affected soils because of its saline and environmental tolerance [53], but it sometimes causes negative ecological impacts for rapidly invading lands [52]. The Leucaena collinsii Britton and Rose species (known as guaje) is a deciduous tree of the Fabaceae family, distributed from southern Mexico to Guatemala and is widely used as forage in agroforestry systems [20].

Evaluation of Survivorship and Vegetative Growth of Species in Field Conditions
The survival rate, height, diameter, and cover of the species were evaluated in field conditions to obtain a quantitative measure of their adaptation and vegetative growth under the influence of the edaphoclimatic factors of the site.

Survival Curve Estimation: The Kaplan-Meier Estimator
The survival rate was estimated with the Kaplan-Meier estimator, , also known as the limit product estimator, which is a nonparametric statistical method in which censoring is considered when studying the period that is necessary for an event to occur. This estimator is used when the individual times of the censored and the uncensored are

Target Species
The Atriplex canescens (Pursh) Nutt. species known as Chamizo (fourwing saltbush) is a shrub in the Amaranthaceae family native to the arid zones of North America, distributed from Canada to central Mexico [50]. This species is an important food source, both for both domestic animals (mainly cows and goats) and for wild animals [51], and it is used for reforestation, soil conservation, and soil recovery [21]. Cynodon dactylon (L.) Pers. (Bermuda grass) is a perennial grass in the Poaceae family native to tropical Africa, Eurasia, India, and Malaysia and is widely distributed in tropical and subtropical regions in the world [52]. This species is used as a source of food for livestock as well as a phytoremediation tool for salt-affected soils because of its saline and environmental tolerance [53], but it sometimes causes negative ecological impacts for rapidly invading lands [52]. The Leucaena collinsii Britton and Rose species (known as guaje) is a deciduous tree of the Fabaceae family, distributed from southern Mexico to Guatemala and is widely used as forage in agroforestry systems [20].

Evaluation of Survivorship and Vegetative Growth of Species in Field Conditions
The survival rate, height, diameter, and cover of the species were evaluated in field conditions to obtain a quantitative measure of their adaptation and vegetative growth under the influence of the edaphoclimatic factors of the site.

Survival Curve Estimation: The Kaplan-Meier Estimator
The survival rate was estimated with the Kaplan-Meier estimator,Ŝ KM , also known as the limit product estimator, which is a nonparametric statistical method in which censoring is considered when studying the period that is necessary for an event to occur. This estimator is used when the individual times of the censored and the uncensored are known, thereby calculating survival rates whenever an individual dies or reaches their follow-up date.Ŝ KM = ∏ t i <t , where r(t i ) is the live plants, d(t i ) is the dead plants and t i < t is the time at which the measurement was made [42,54,55]. In accordance with Kaplan and Meier [54], the survival function is defined as follows: where S(t) is the probability of a death occurring at a time T at least equal to the time t [54].
To this end, the status of each plant (alive or dead) at the end of the evaluation period was considered, as well as the length of its life (months). Considering F(t) to be the distribution function of T, then S(t) = 1 − F(t), indicating that S(t) is decreasing. The rate at which S(t) decreases at a particular instant indicates the probability (risk) that the individual will die at that moment [54]. The analysis was conducted using SPSS v25 software.

Quantification of Growth Analysis: Plant Height, Stem Diameter, and Canopy Cover
The species height (cm) was measured from the base of the stem to the apex of the dominant vertical shoot, whereas the diameter of the stem (cm) was measured at ground level. The canopy cover (%) was calculated using the quadrat method. The initial height (cm), the diameter at ground level (mm), and the canopy cover were obtained one month after out-planting. The relative annual growth rate, termed as diameter growth rate (DGR, cm month −1 ) for the diameter at ground level, and the height growth rate (HGR, cm month −1 ) for height and plant cover (CGR, % month −1 ) for cover were calculated following the equation [56]: where W represents either diameter, height, or cover at t 1 at the beginning of the experiment and t 2 at the end of the experiment.

Species Presence Records
Records of the presence of the species A. canescens, C. dactylon, and L. collinsii were obtained from the databases of the National Biodiversity Information System (SNIB) and verified with the databases of the Global Biodiversity Information Facility (GBIF) and the National Herbarium of Mexico (MEXU). The spatial records for each species were subjected to a randomness test using R [57], ILWIS 3.43 (http://52north.org/ilwis, accessed on 16 May 2019) and QGIS 3.12 (https://qgis.org/es/site/, accessed on 16 May 2019), in order to reduce or eliminate their inherent autocorrelation. The test is used to determine whether a sample or data set follows a pattern or is random [58]. If the randomness test is positive, 75% of the records are used for training and 25% for validating the model; otherwise, a pattern analysis is applied by estimating the distance at which it is possible to find one species record with the greatest probability [59,60]. Pattern analysis was used to determine the distance at which the records do not show spatial autocorrelation [61], using ILWIS 3.43 software (http://52north.org/ilwis, accessed on 16 May 2019).

Edaphoclimatic Variables
To calculate the planting potential of the three species, 19 edaphoclimatic variables were extracted for each point of occurrence, homogenized in raster format with 30 arcsecond resolutions. Five climatic variables were obtained from the WorldClim Version 1 platform (http://www.worldclim.org/, accessed on 8 December 2019) [62]: annual mean temperature (BIO1), isothermality (BIO3), seasonal temperature (BIO4), annual precipitation (BIO12), and seasonal precipitation (BIO15). Three evapotranspiration variables generated with the Turc model [63]: real annual evapotranspiration (ETR A ), real evapotranspiration of the humid months (ETR H ), and real evapotranspiration of the dry months (ETR S ). Evapotranspiration was calculated using the Turc's model [63]: where P = total annual precipitation (mm), L = 300 + 25T + 0.05T 3 , and T = average annual temperature ( • C). The other used climatic variables were two seasonal variables corresponding to the humid (from May to October) and dry (from November to April) months: precipitation of the humid months (PP H ), precipitation of the dry months (PP S ), mean temperature of the humid months (T H ), and mean temperature of the dry months (T S ). The three edaphic variables related to soil properties used in the model were: organic matter (MO), electrical conductivity (CE), and hydrogen potential (pH), generated by the Soil Profiles 1:250,000 Series II data set [59]. In addition, there were used two topographic variables generated from a Digital Elevation Model, with a pixel size of 120 m: terrain orientation (ASP) and elevation (ELEV), and two normalized vegetation indices obtained from remote sensing data (MODIS): normalized vegetation index for the most humid month (INV H August 2019) and normalized vegetation index for the driest month (INVs January 2019).

Adjustment and Evaluation of the Predictive Model
The linear correlation (r) between the edaphoclimatic variables was verified using the corrplot package in the R software [57]. Variables with correlations greater than or equal to 0.80 were eliminated to avoid multicollinearity in the model. Based on the occurrence points and selected edaphoclimatic variables, the species aptitude probabilities were estimated with MaxEnt, version 4.4.1, which is a machine learning method developed for species distribution modelling and for mapping their habitat suitability with remarkably high predictive accuracy [29,64,65]; it was chosen for their good performance, use, and operation [46]. The estimated probabilities of species occurrence range from 0 to 1. The configuration in the MaxEnt algorithm was set by default [66]. The random seed option was selected so that each repetition would use a different set of test points. The model was submitted to 10 repetitions, the subsampling was performed for each repetition, and the 75 percent of the data set was used for training and 25 percent for testing. The convergence threshold was set at 0.00001, using 500 interactions and 10,000 background points, and the results were based on an average of 10 repetitions. The predictive capability of the models was assessed by the area under the curve (AUC) of the receiver operating characteristics (ROC) that measures the ability of the model to discriminate between sites where a species is present versus sites where it is absent. AUC values greater than 0.9 are considered models with high predictive capability, between 0.6 and 0.9 with good predictive capability and less than 0.5 predictive capability lower than expected at randomness [67,68]. The contribution of each edaphoclimatic variable to the model's area under the receiver operator curve (AUC) was quantified using the Jackknife test, which involves removing the variables using a step-by-step procedure until a minimal model is left, consisting only of variables that increase model performance [69].
The MaxEnt model results were imported into the ArcGIS software. The final maps for the species planting potential comprised values between 0 and 1 of habitat suitability. According to Qin et al. [27] and Yang et al. [70], values greater than 0.6 are considered sites with high planting potential, values between 0.4 and 0.6 sites with good planting potential, and values lower than 0.4 sites with low planting potential.

Survival Curve Estimation: The Kaplan-Meier Method
The results indicate that 90.7% of the shrubs and the grass established in the experimental plot remained alive after being evaluated for 15 and 3 months, respectively ( Table 1). As a result of censored data, which are plants that did not succumb to the event death during the analysis, L. collinsii had the highest estimates of surviving with 100% (Figure 3), whereas A. canescens had the highest estimate of survival risk with 20% ( Figure 4).

Survival Curve Estimation: The Kaplan-Meier Method
The results indicate that 90.7% of the shrubs and the grass established in the experimental plot remained alive after being evaluated for 15 and 3 months, respectively ( Table  1). As a result of censored data, which are plants that did not succumb to the event death during the analysis, L. collinsii had the highest estimates of surviving with 100% (Figure 3), whereas A. canescens had the highest estimate of survival risk with 20% ( Figure 4).   In Table 2, it is showed the average growth rates of plant height, stem diameter, and canopy cover achieved by each of the species during the evaluation period. A. canescens and L. collinsii species exhibited growth rates of 0.71 and 1.78 cm month −1 in height, 0.04 and 0.11 cm month −1 in diameter, and 0.54 and 0.15% month −1 in cover, respectively. C. dactylon obtained average rates of growth of height and cover of 0.94 cm month −1 and 7.47% month −1 , respectively, after transplanting.

Quantification of Growth Analysis: Plant Height, Stem Diameter, and Canopy Cover
In Table 2, it is showed the average growth rates of plant height, stem diameter, and canopy cover achieved by each of the species during the evaluation period. A. canescens and L. collinsii species exhibited growth rates of 0.71 and 1.78 cm month −1 in height, 0.04 and 0.11 cm month −1 in diameter, and 0.54 and 0.15% month −1 in cover, respectively. C. dactylon obtained average rates of growth of height and cover of 0.94 cm month −1 and 7.47% month −1 , respectively, after transplanting.

Species Presence Records
Four hundred and twenty-eight records were found for A. canescens, 1452 for C. dactylon, and 74 for L. collinsii. The records of the species A. canescens were mostly located in the states of Coahuila, San Luis Potosi, Baja California, Sonora, Chihuahua, Baja California Sur, Nuevo León, Zacatecas, and Mexico State. The records of C. dactylon were found in all the states of Mexico, and those of L. collinsii were found in the states of Chiapas, Mexico, and Oaxaca. The results of the randomization test applied to the records for each species indicate that the records for the three species show aggregate spatial. The records of C. dactylon and L. collinsii show moderate aggregation at distances less than 30 km on average (≈0.28 degrees), being the probability of finding a record high (G (r) ≈ 0.9). On According to the pattern analysis, the number of final records required to avoid overfitting in the models was 150 for A. canescens, 486 for C. dactylon, and 36 for L. collinsii, which when applied, satisfied their randomness test.

Adjustment and Evaluation of the Predictive Model
Based on the ROC analysis, the AUC values for evaluating the predictive capability of the models for A. canescens and L. collinsii were 0.845 and 0.992, respectively, for training data, and 0.867 and 0.997, respectively, for test data. The AUC value for evaluating the predictive capability of the model for C. dactylon was 0.721 for the training data and 0.697 for the test data. As a result of correlation analysis and variance inflation factor analyses, which demonstrate the degree of multicollinearity in a set of multiple regression variables [71],  Figures 5 and 6 show the optimum planting areas for A. canescens and C. dactylon species in the Tula watershed, which were good (>0.4) and covered 355 km 2 (34%) and 477 km 2 (46%) of the total watershed area, respectively. The optimal planting areas for L. collinsii were low (<0.4), and they covered 100 percent of the total watershed (1037.66 km 2 , Figure 7). The semiarid temperate regions (72.96%) hold more potential for A. Canescens planting than the humid areas (20.53%) of the watershed, whereas the optimal areas for planting C. dactylon were equally spread throughout the temperate semiarid regions (46.20%) and subhumid climate zones (46%). The species L. collinsi did not present optimal planting sites within the Tula watershed. km 2 (46%) of the total watershed area, respectively. The optimal planting areas for L. collinsii were low (<0.4), and they covered 100 percent of the total watershed (1037.66 km 2 , Figure 7). The semiarid temperate regions (72.96%) hold more potential for A. Canescens planting than the humid areas (20.53%) of the watershed, whereas the optimal areas for planting C. dactylon were equally spread throughout the temperate semiarid regions (46.20%) and subhumid climate zones (46%). The species L. collinsi did not present optimal planting sites within the Tula watershed.

Evaluation of Survivorship and Vegetative Growth of Species in Field Conditions
The average survival rate for the three species was excellent (90.7%) at the experimental site. However, Omary [72] points out that plant survival can be affected by soil properties such as moisture, temperature, pH, electrical conductivity, and nutrient content as well as the planting or transplanting method [73]. In this study, the experimental site, where the three species were transplanted, has the same soil as substrate, thus the differences in their survival rates were probably related to environmental conditions and topographic factors such as slope and elevation [74].
The survival rate for the species A. Canescens was acceptable (80%) a year after transplanting. This value is as the report of Saucedo and Chacón [75] who estimated a survival rate of 71 to 85% under dry weather conditions in Chihuahua for the same evaluation period. However, this percentage was higher than the survival rate (65.4% and 63.9%) reported by Rios-Saucedo et al. [76] for degraded areas of Chihuahua, Mexico, 60 and 150 days following planting, both under a silvopastoral system. For the species A. canescens, Saucedo and Chacón [75] reported an average height of 20.3 cm and 49.9 cm after six and fifteen months of growth, respectively, under planting frames of 2 m around each plant and dry climatic conditions in Chihuahua, Mexico. In this study, the average height obtained was 23.8 cm after 12 months of planting with a HGR of 0.71 cm month −1 . This height is lower than the one found by Saucedo and Chacón [75] fifteen months after planting but slightly higher than that obtained six months after planting, indicating that the vegetative growth of A. canescens is slower under the soil and climatic conditions of the study area. These differences in height could have been caused by the low temperatures of the winter months, as well as by the lack of moisture in the soil, both of which negatively impact the growth of roots in young A. cenescens (two years old), whose roots are underdeveloped [75,77].
The species C. dactylon showed excellent survival rate (92%) after three months of planting. This value was higher than the 75 percent reported by Laurencena et al. [78] for tropical and temperate regions, with an annual mean temperature of 18.1 • C, an annual precipitation of 947.6 mm per year and a molisol soil with a superficial horizon of silt loam texture. The average height of C. dactylon species was 3.3 cm for a three-month evaluation period with a HGR of 0.94 cm month −1 . According to Hernandez et al. [79], Bermuda grass Tifton 68, a fast-growing variety of C. dactylon, requires three months to reach 35 cm in height before it can be grazed or cut (establishment period). Thus, this study indicated that the growth rate of C. dactylon was influenced by the genotype, its origin, and its adaptation to soil and environmental conditions. The survival rate for the L. collinsii species was excellent (100%) after 12 months of transplanting. Osuna-Ceja et al. [80] reported similar results for another species within the genus, L. Leucocephala, in an agroforestry system with cactus plants as living barriers, in a semi-arid highland valley of north-central Mexico. This valley, located at an elevation of 2049 m asl, has an average precipitation of 344 mm during the growth cycle of the crop, an average temperature of 16.3 • C, and a soil depth of less than 45 cm. Within its arable layer, the soil has a sandy loam texture with less than 1% organic matter, a slope of 2%, and a pH of 6.6. Membreño [81] also reported similar survival rates (94%) with an average height of 169 cm and a stem diameter of 2.07 cm for the variety L. collinsii 45-85 established in a transition zone between tropical dry forest and tropical subhumid forest in Nicaragua. The average height and diameter in this study for the species L. collinsii were 52.4 cm and 3.8 cm, with a HGR and DGR of 1.78 cm month −1 and 0.11 cm month −1 respectively, twelve months after planting. The HGR and CGR in this study increased over time, but they were low when compared to those in other studies. According to Burner et al. [82], these results are the consequence of a lack of moisture, which causes a reduction in branch elongation, plant growth rate, and net accumulation rate [83]. Although the L. collinsii species reached a lower height and a larger stem diameter compared to other studies, its capacity to adapt to the site was high due to its pivotal root system that can reach nutrients and moisture from deeper soil layers causing a high regrowth capacity in times of lack of moisture. Although the L. collinsii species achieved a slightly lower plant height and stem diameter compared to other studies, the adaptability of the species to the site was high due to its pivotal root system that can reach nutrients and moisture from deeper soil layers, thereby supporting a high growth rate even in the absence of moisture [84]; probably due to digging a hole and adding substrate. Its adaptation was likewise due to its average foliar retention (40-60%), a factor that is relevant in environments with dry periods of six months or longer [85].

Determination of the Potential Planting Areas for the Species
The models were generally accurate in predicting the planting areas of the three species, as their training and testing AUC values were all greater than or equal to 0.70, indicating that they can differentiate between sites with high and low planting potential [68,69]. In addition, the models indicate that climatically more suitable areas or areas with high planting potential are positively correlated with the environmental conditions found in the natural habitats of these species [86]. The prediction capability of the C. dactylon model was low (AUC = 0.7) in comparison with those of the L. collinsi and A. canescens models (AUC = 0.9), which is explained by their wide geographic distribution [70,87].
The species A. canescens presented great potential for planting in the semiarid climate regions of the study area (>0.4). These regions are characterized by average temperatures of 14 • C and an average precipitation of 500-600 mm per year. The predominant soils are haplic phaeozem (high organic matter, shallow, stony, and extremely unstable) and eutrophic cambisol (extremely productive soils with a high degree of susceptibility to erosion), both of medium texture (loamy) [19,[44][45][46]. The high planting potential is strongly related to precipitation, soil properties and the normalized vegetation index for the driest month (specifically: PP H , 52.9%, soil pH, 12.5% and INV S , 9.3%). Normalized vegetation index and climatic variables, especially precipitation, must be examined with respect to the seasonal variations in vegetation, as water availability changes considerably during the various phases of vegetation development [88]. Wang et al. [88] demonstrated that INV has a strong linear relationship with precipitation during the growing season. However, the INV did not have a significant relationship with January precipitation due to insufficient soil moisture because of high temperatures, evapotranspiration, and drought. The high potential of planting in the region is due primarily to the seasonality of the precipitation and the high evaporation which allows the plants to achieve optimal growth [21,89].
The species C. dactylon showed great planting potential (>0.40), with suitable planting sites evenly distributed throughout the Tula watershed, influenced strongly by temperature, elevation, and precipitation (specifically: T S , 15.4%, IVN S , 15%, ELEV, 12%, BIO15, 11.9%, BIO01, 10.2%, y INV H , 7.8%). The high suitability pattern related to climate between sites is particularly important for the planting and vegetative growth of C. dactylon as warm temperatures are among the most important factors that promote its adequate development [78]. Although the C. dactylon species is generally found in low-elevation areas, the average altitude (2778 m asl) in this watershed was favorable for its growth, corroborating the fact that it is a cosmopolitan plant [52].
The species L. collinsii did not have potential to be planted (<0.4) in either of the two climatic zones, temperate semi-dry or sub-humid temperate, within the watershed. The temperatures in these areas are typically around 14 • C with precipitation of 500-600 mm and 600 to 1000 mm per year on average. The predominant soils are the haplic phaeozem and the eutrophic cambisol, both of medium texture (loamy), but due to their management they are highly eroded and shallow soils, as well as the pelic vertisol, a fertile soil with high shrink-swell potential with cracks developing under dry weather. Planting potential is strongly affected by seasonal precipitation, elevation, and seasonal temperature (specifically: PP H , 32.7%, ELEV, 11.1%, BIO15, 10.6%, T S , 10.1%, and BIO4, 9.4%). This pattern of suitability for the species L. collinsi in relation to the climatic environment of the watershed is particularly important for its establishment and vegetative development, since two major factors constraining its expansion are low precipitation and elevation [20,90,91], despite its high tolerance to drought [92]. The mean elevation of the watershed is 2878 m above sea level, which is considerably higher than the range of elevations between 400 and 900 m where the species can develop successfully [20].
Low precipitation also influences the availability of nutrients [93] that alter the differentiation of new tissues (regrowth meristems) and the growth of those already in existence, as well as a decrease in foliage area, resulting in a faster rate of leaf senescence [94,95]. The Jackknife test confirmed that the effect of precipitation and temperature variables, more specifically their seasonality in the plantation modeling, was significant in the three species. Recent studies have shown that seasonal temperature variations can be used to predict suitable planting sites areas [86,96,97] because photosynthetic activity is often influenced by them more than by extreme temperatures [98].

Conclusions
Environmental modeling using the Kaplan-Meier estimator and the MaxEnt algorithm made it possible to determine the survivorship within potential plantation areas of the three non-native species within the Tula watershed; the simulation models can be used to plan the introduction of non-native species. The plant height, stem diameter, canopy cove, and survivorship were affected by environmental factors, mainly related with precipitation, temperature, and soil fertility, resulting in that the morphological characteristics of the plants were not fully expressed. A. canescens shrubs and C. dactylon grass presented areas with great planting potential and excellent survival rates, respectively, suggesting that the Tula watershed offers an optimum environment for their plant growth and development. Despite having an excellent survival rate in the experimental plot, the potential distribution model for L. collinsii suggests that there are no areas with optimal environmental conditions for its growth. The results of this study provide a solid basis for selecting species by evaluating their vegetative growth and survivorship, as well as their potential planting sites, to recommend agroforestry as a land management system on degraded soils with species that