Phenotypic, Geological, and Climatic Spatio-Temporal Analyses of an Exotic Grevillea robusta in the Northwestern Himalayas

Abstract

The last five decades (since 1980) have witnessed the introduction of exotic trees as a popular practice in India to fulfill the demand of forest-based products for utilization in afforestation programmes. This study examines the distribution and habitat suitability of exotic Grevillea robusta trees in the northwestern Himalayas (state: Uttarakhand), focusing on the interaction between G. robusta and abiotic factors, such as climate, soil, and habitat suitability. This multipurpose agroforestry species is mainly grown by farmers as a boundary tree, windbreak, or shelterbelt and among intercrops on small farms in agroforestry systems worldwide. The results indicate that phenotypic plasticity is determined by tree height and diameter, indicating a higher frequency of young and adult trees. The study also highlights spatio-temporal modeling coupled with geological analysis to address the current distribution pattern and future habitat suitability range through MaxEnt modeling. The AUC ranged from 0.793 ± 3.6 (RCP 6.0_70) to 0.836 ± 0.008 (current) with statistical measures, such as K (0.216), NMI (0.240), and TSS (0.686), revealing the high accuracy of the model output. The variables, which include the minimum temperature of the coldest month (Bio 6), the slope (Slo), the mean temperature of the driest quarter (Bio 9), and the precipitation of the driest quarter (Bio 17), contribute significantly to the prediction of the distribution of the species in the Himalayan state. The model predicts a significant habitat suitability range for G. robusta based on bio-climatic variables, covering an area of approximately ~1641 km² with maximal occurrence in Pauri (~321 km²) and Almora (~317 km²). Notably, the future prediction scenario corroborates with the regions of Tons (Upper Yamuna, Uttarkashi), Kalsi (Mussoorie, Dehradun), the Kedarnath Wildlife Sanctuary, and the Badrinath Forest Division for the potentially suitable areas. The climate was found to have a strong influence on the species’ distribution, as evidenced by its correlation with the Köppen–Geiger climate classification (KGCC) map. While the species demonstrated adaptability, its occurrence showed a high correlation with bedrocks containing an elevated iron content. Furthermore, the study also provides the first trees outside forests (TOF) map of G. robusta in the region, as well as insight into its future habitat suitability.

1. Introduction

Worldwide, forests and agroforestry areas are undergoing consistent plantation of exotic species to fulfill various economic and environmental demands [1,2]. Grevillea robusta A. Cunn. Ex R. Br. (silver oak; family: Proteaceae), for example, is a species whose wood has economic and ecological potential [3,4,5]. Economically, it is grown as a shade tree for tea and coffee plantations, timber, firewood, and poles, with the leaves being used for mulch. Ecologically, the species also stabilizes soils and enhances water and nutrient uptake. The roots help it to grow in low-fertility soils [6] and recycle nutrients from the deep soil surface to the upper soil surface through the decomposition of leaf litter [7].

The species is endemic to Queensland, i.e., eastern Australia, and is widely distributed in southeast Asia, Africa, and the southern part of North America and South America [8]. In the last three decades, G. robusta was introduced from Australia to all regions of the world, viz., China, Eritrea, Ethiopia, India, Indonesia, Jamaica, Kenya, Laos, Malawi, Malaysia, Mauritius, Nepal, Pakistan, the Philippines, South Africa, Sri Lanka, Tanzania, Uganda, the United States of America, Vietnam, Zambia, and Zimbabwe [9,10,11,12,13,14,15,16]. In India, the species is distributed in approx. more than 20 states of India, such as Andhra Pradesh, Arunachal Pradesh, Assam, Delhi, Haryana, Himachal Pradesh, Kerala, Karnataka, Madhya Pradesh, Maharashtra, Nagaland, Orissa, Punjab, Tamil Nadu, Telangana, Uttarakhand, Uttar Pradesh, and West Bengal (World Agroforestry Centre 2002), showing its reproductive fitness across the varied agroclimatic zones of the country.

In an uneven landscape like the Himalayas, basic information on the pattern of species distribution and a genetic resource assessment for trees outside forests (TOFs) map is lacking and needs sound statistical assessment. In 2000, the United Nations Economic Commission for Europe (UNECE) and the Food and Agriculture Organization (FAO) TOFs were defined as “trees on land that meets the definitions of forest and of other wooded lands except that the area is less than 0.5 ha and the width is less than 20 m”. Importantly, there is a need to understand and predict the relative habitat suitability range for any exotic species before their introduction [17,18]. A large number of statistical models are co-currently in use to simulate the spatial distribution of plant species [19,20], the habitat suitability of the species [21,22], and the spatial pattern of species diversity [23], as well as to predict the impact of climate change [24,25]. Grevillea robusta has a fast-growing nature and an acclimatization ability across the varied environmental conditions, which reveals the species’ adaptability and suitability to the vast distributional range [3,5]. Importantly, life-history traits (such as fast growth, high propagule production, hybridization, etc.) and favorable environmental conditions (such as low and intermediate shade tolerance, site and biotic factors, canopy disturbances, etc.) are the main factors controlling the distribution range of exotic plants [1].

In recent years, policymakers have become more conscious of the importance of TOFs in terms of the possible economic roles and contributions to human wellbeing through ecosystem services. As a result, tree-based resources are frequently evaluated in forest monitoring using a remote sensing and geographic information system (RS&GIS) [26]. For instance, the TOFs’ biomass was demonstrated to play a substantial role in regional and national bio-economies when analyzing TOFs data across a 13-county area of Michigan (USA) [27]. As the world’s population grows, TOFs are destined to become a major component of natural resources, require mapping to determine the expanse and extent of the distribution and monitoring to improve green cover [28,29,30].

Importantly, the assessment and monitoring of various TOFs resources necessitate the accuracy of data and variables [FCM, forest type map (FTM), non-forest areas, etc.]. High-resolution data were used in conjunction with an object-based image analysis (OBIA) technique to meet the changing demands of resource assessment and monitoring. Satellite imaging, such as the LANDSAT series, SENTINEL, QuickBird, RapidEye, IKONOS, and Google Earth, ranges in resolution from low to high by the European Space Agency, the German Aerospace Center, and the National Aeronautics and Space Administration. Satellites such as CARTOSAT-I, Indian Remote Sensing (IRS-P6), linear imaging self-scanning sensor (LISS-III), LISS-IV, and others, demonstrate simulated data at various spatial resolutions with more consistent and reliable information on TOFs by the Indian Space Research Organization [31,32,33,34,35]. For the classification of linear patches of TOFs along roadside and field bunds, the IKONOS fused image (MSS and PAN) showed itself to be superior [36]. These effectively provide credible information on urban greens and forest cover in environments that have been heavily influenced by humans [37,38,39,40,41]. The RS&GIS detects changes in forest cover and TOFs concerning land use/land cover (LULC) primarily using historical and recent satellite photos (LANDSAT and RapidEye), which are supplemented by ground verification [42]. As a result, this method appears to be practical for identifying and stratifying TOFs.

Woodlands are usually exposed to heterogeneity in the environment due to various factors (climate change, edaphic, topographical and land use changes), and differences are observed in terms of the phenotypic response, as shown by the morphological traits, viz. height, diameter, branching pattern, abiotic and biotic interference, which reflect the varied adaptive response. In silviculture and forestry practice, it is necessary to predict the variation in morphological traits [43] as the size class distribution reveals information about the stand structures (SS) and its dynamics. This type of study reveals the demographic stability and phenotypic plasticity of a species in an area. Phenotypic plasticity is considered one of the major means by which exotics, such as Eucalyptus sp., Populus sp., and Grevillea sp., cope with the variation of factors specific to each locality. Individuals of a species often experience different environmental conditions, which often require adaptation in a particular habitat to maximize fitness [44]. Notably, a genetic resources assessment includes the relative abundance, richness, size variation, and spatial variation, which might facilitate reliable, precise, and verifiable evidence on the expected size class of trees within its ecological niche [45,46], commonly described by frequency distribution and expressed on a unit area basis [47]. According to Pastorino et al. (2021) [48], “Genetic resources from any biological organism-complete genomes, genes, or even portions of genes that have actual or potential relevance to humans”. Thus, in tree species, it becomes necessary to categorize and characterize the morphological traits in which phenotypic plasticity may play a significant role to predict global climate change and the ecological consequences. For instance, land use change and invasiveness at an ecosystem level can be assessed. As individual genotypes have the potential for adaptive variance, plasticity in plants possibly broadens ecological ranges and reduces the impact of selection [49,50].

As per the insufficient literature available in the Indian context, G. robusta has not yet been studied for habitat suitability, phenotypic plasticity, and quantification of genetic resources, particularly in the Himalayas. It is hypothesized that the generation of baseline data and information through species distribution modeling (SDM) is an important tool to assess the potential habitat suitability range and key factors determining the phenotypic plasticity of this species; therefore, this work represents the characterization of G. robusta in the northwestern Himalayas by targeting the following objectives: (i) quantification of genetic resources through size class distribution and assessment of the phenotypic plasticity on the basis of morphological parameters; (ii) exploration and mapping of the current extent of the distribution range using a MaxEnt modeling approach and plotting on the Köppen–Geiger climate classification (KGCC) and geological maps; and (iii) prediction of the spatio-temporal habitat suitability range of G. robusta in accordance with future climate change scenarios. Holistically, the study provides an opportunity for the suitable implication to map a heterogenous formation of TOFs, and direction for exotic tree management in the landscape ecosystem of the Himalayas.

2. Materials and Methods

2.1. Field Survey and Sampling Strategies

Grevillea robusta was planted in an inhabited site, most likely by forest officials and human interventions; therefore, the field surveys were conducted throughout the road networks (rural, state, and national highways), government institutions, schools, hospitals, industrial areas, etc. (planted as avenue trees), in the 13 districts representing the state of Uttarakhand of the northwestern Himalayas during 2018–2022 (Figure 1).

The study area lies between N 28°43′31.288″ to 31°27′07.908″ latitude and E 77°34′26.805″ to 81°02′23.228″ longitude with an elevation ranging from 150 to 1934 m msl (Figure 2).

A multiphase random sampling methodology was applied. Initially, extensive surveys were conducted to determine the species occurrence areas. The sampling followed linear transects (quadrats) of 100 m laid out, where at least one tree was present, and a width of 10 m was considered for each transect forest survey of India (FSI 2010). In the second phase, geospatial data (latitude, longitude, altitude) was recorded with a GARMIN Etrex 10 handheld global positioning system (GPS) device in the degree–minute–decimal format with a positional accuracy of 5–8 m for each sampling plot. The trees were counted all along the roadsides, institutional lands, and parks. This unbiased sampling approach in different stages was adopted to reduce the cost, but it also allowed us to focus on the in-field operations around the selected points for precise estimates. The methodology also provides the potentially accurate estimate of species distributions with minimum errors [51]. A total of ~1557 quadrats were recorded from all thirteen districts of Uttarakhand.

2.2. Size Class Distribution

The data on ~13,000 trees representing ~1557 quadrats were recorded by assuming that all the trees were present in the same plot with same size characteristics (as the observed trees were planted ones). The morphological parameters, such as tree height (TH, m) and clear bole height (CBH, m) measured by a Bluemish altimeter, the girth at breast height (GBH, m) and crown width (CW, m) measured by meter tape, and branch angle [BA, (°)] measured by a visual superimposition of a protractor were also recorded. Furthermore, the diameter at breast height (DBH, m) was calculated from GBH by dividing by 3.14 (G = πD, where G = girth, D = diameter, π = 3.14).

During recording of the field data, biotic factors, such as the incidence and severity of diseases (DI) or insect or pest (IPI) attacks on the trees were also taken into consideration. These symptoms were analyzed by ocular observation and noted down by a given scalar method to determine the actual severity (

Table 1. Scoring pattern of individual tree of G. robusta for in-field observation.

Sl. No.	Category		Score
	Disease Severity	Insect/Pests Incidences
1.	No infection	No incidence	1
2.	1–15%	1–15%	2
3.	15–30%	15–30%	3
4.	30–45%	30–45%	4
5.	>45%	>45%	5

).

Each tree was given a score according to the extent of DI and IPI, and the final severity was evaluated as per the standard formula [52,53] (Koyshibayev and Muminjanov 2016). The overall size class distribution of G. robusta for variables TH and DBH were classified by dividing the samples into four size classes in accordance with their altitudinal gradient, i.e., <601 m, 601–1050 m, 1051–1500 m, and >1500 m msl, and were then classified into four frequency classes representing juvenile, young, adult, and mature trees on the basis of DBH [54]. Similar classes were determined for the other morphological parameters, such as CBH, CW, BA, DI, and IPI.

2.3. Climatic and Geological Mapping for Current Distribution and Future Prediction

The prediction mapping scheme for G. robusta is shown in Supplementary Figure S1. The forest cover map of the study area was generated by using SENTINEL datasets of Uttarakhand. The output was generated into two major classes (forest and non-forest), but we used the sample data of the non-forest classes map for exotic G. robusta. For this, Google Earth™ imagery was used to vectorize the road network, and the vector product was saved in a Keyhole Markup Language (KML) file format. Furthermore, the KML file was converted into a shapefile using ArcGIS Ver. 9.3 software.

Additionally, the bioclimatic variables were used in SDM with an average monthly climate data for the minimum, maximum, and mean temperature and precipitation over 1970–2000, with a 30 s (~1 km²) resolution of WorldClim Ver. 2 [55]. To create the layers of elevation, slope, and aspect, a scene from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) sensor and a global digital elevation model (GDEM) were used [56]. Direct normal irradiance (DNI) provides the monthly average, annual average, and daily total solar resource averaged over the surface cells of 0.038 degrees in both latitude and longitude, or nominally 4 km in size. Further, the direct normal irradiance (DNI) was downloaded from the Solar Energy Centre, Ministry of New and Renewable Energy (MNRE), Government of India (GoI), New Delhi (https://maps.nrel.gov/nsrdb-viewer/; 25 October 2022).

Finally, maximum entropy (MaxEnt) modeling was used to estimate the current distribution and identify the potential future habitat suitability areas in the Uttarakhand using the MaxEnt Ver. 3.1 program [57]. The geospatial parameters and environmental variables, i.e., iso-thermality (Bio 3), min temperature of the coldest month (Bio 6), mean temperature of driest quarter (Bio 9), precipitation of driest month (Bio 14), precipitation seasonality (Bio 15), precipitation of driest quarter (Bio 17), slope, and aspect, have maximum permutational contribution in prediction and were used in 100 replicates. Important model features, such as linear (0.050), quadratic (0.050), hinge (0.500), categorical (0.250), and threshold (1.000), were used with 10,000 background points. To decrease model over-prediction or overfitting, the regularization multiplier (RM) value was set to 0.1 with 5000 iterations and the rest of the settings were kept as default [58,59,60,61].

The outcome of the model depends on the area under receiver operating characteristic (ROC) curve (AUC) and the jackknife test. The AUC determines how the model runs and its prediction significance, whereas the jackknife test determines the influence and contribution of each variable used in prediction modeling [57,62,63]. The current prediction distribution had a probability range from 0 to 1, which were regrouped, and those with a value >0.7 were used in the present study. The distribution map area was overlaid on the KGCC (1976–2000) derived from Kottek et al. (2006), [64] and the geological map (geological survey of India; http://bhukosh.gsi.gov.in; 15 November 2022). The objective of comparing the current distribution to the geological map of Uttarakhand is to investigate the possible influence and correlation of bedrock types in the habitat suitability range of G. robusta. With respect to the scale of this study, a map of the current (Figure 2) and future maximum possible habitat suitability occurrence prediction of RCP 8.5_2050 (1000 m buffer) was overlaid on a simplified geological map (1:2M scale, geological survey of India) to avoid the geological complexity that is typical of orogenic regions like the Himalayas.

Effectively, the MaxEnt was further used for future prediction, where the data of two time periods, i.e., 2041–2060 (2050s) and 2061–2080 (2070s), signify four representative concentration pathways (RCPs), i.e., 2.5, 4.5, 6.0, and 8.5 were downloaded from WorldClim (http://www.worldclim.org; 18 February 2022). The RCPs data were used with the same set of eight bioclimatic variables, and the output > 0.7 was masked by a 500 m and 1000 m buffer to generate the probable habitat suitability prediction map in accordance with the climate change scenarios. Additionally, the accuracy of the habitat suitability map (generated through the MaxEnt model output) was evaluated by estimating the kappa coefficient (K), normalized mutual information (NMI) n(s), true skill statistics (TSS), etc., using the derivatives of the confusion matrix [65,66].

3. Results

3.1. Geographical Distribution and Stand Structure (SS)

The field data revealed that G. robusta distribution ranges from 165 to 1957 m amsl, and occurs mainly in roadside plantations, community park land, schools, and official premises. A widespread distributional range with respect to altitude and latitude was also observed at the landslide-prone areas of national and state highways (district Uttarkashi and Chamoli), where the species has reclaimed and stabilized the slope. The tree is majorly found in the planted strips and is associated with other roadside planted species, such as Alstonia scholaris, Cassia fistula, Holoptelia integrefolia, Saraca asoka, and in some places with Pinus roxburghii.

In total, ~13,000 trees were counted for morphological parameters, and each tree was distributed into one of four categories (juvenile, young, adult, and mature trees) along the altitudinal gradient (<601–>1500 m msl). For the stand structure (SS) size class distribution analysis, parameters, viz., TH and DBH (derived from GBH; Figure 3b), were used. For all the altitudinal classes, the frequency distribution of the height class (ranging from 9.1–27.0 m) and the diameter class (ranged from 0.24–0.44 m) revealed that young trees and adult trees were the most frequent (Figure 3a–c). Other parameters, such as CBH, CW, and BA, were also measured to determine the characteristic superiority and potential to adapt in varied localities. The CBH class revealed that up to 5 m in tree height, the maximum number of trees had a good, clear, cylindrical bole (Figure 3d). A large number of trees occupies the maximum range (>3.0 m) of CW, which shows G. robusta dominance for height superiority (Figure 3e). The branch angle class indicates that a maximum number of individuals are near to a right angle, i.e., >70.0° (Figure 3f). Notably, the DI and IPI for each tree revealed that there was a low to very low level of biotic severity in the studied stand (Figure 3g–h).

3.2. Current Habitat Suitability Mapping through MaxEnt Modeling

A total of 1167 geo-coordinates were recorded from 13 districts of Uttarakhand to map the geographical distribution (Figure 2). In total, ~821 (70.35%) GPS points were used to test the MaxEnt model, and the remaining were utilized for validation. The performance was evaluated by AUC, and the model with the best performance was used for prediction mapping (

Table 2. The AUC, permutation importance, and the percentage contribution of the different variables.

Time Period, Climate Change Scenario/ AUC Value (Percentage)	Labels	Bio 3	Bio 6	Bio 9	Bio 14	Bio 15	Bio 17	Slop	Asp
	Variables	Iso-Thermality [(Bio 2/ Bio 7) × 100]	Min. Temperature of Coldest Month (Std. Deviation × 100)	Mean Temperature of Driest Quarter	Precipitation of Driest Month	Precipitation Seasonality (Coefficient of Variation)	Precipitation of Driest Quarter	Slope	Aspect
	Scaling Factor	100	100	10	1	1	1
	Units	°C	C of V	°C	mm	mm	mm	0	0
Current/ 83.6 ± 0.8	Percent contribution	2.8	32.3	13.9	2.5	4.5	13.0	15.1	4.8
	Permutation importance	1.8	33.9	14.1	3.9	6.7	18.4	12.7	2.7
RCP 2.6–50/ 82.9 ± 2.8	Percent contribution	6.7	49.4	11.7	3.9	6.6	10.0	9.5	2.2
	Permutation importance	6.0	17.2	48.8	2.8	8.1	10.0	6.4	0.6
RCP 4.5–50/ 82.6 ± 2.6	Percent contribution	5.3	47.6	6.5	1.5	6.1	18.5	12.0	2.4
	Permutation importance	2.2	57.0	16.1	0.9	3.7	14.0	5.6	0.5
RCP 6.0–50/ 82.8 ± 2.5	Percent contribution	8.0	51.5	3.0	2.5	4.4	16.1	11.8	2.7
	Permutation importance	6.3	49.7	5.0	5.1	5.5	19.2	8.2	5.1
RCP 8.5–50/ 82.1 ± 2.7	Percent contribution	5.1	52.8	8.0	1.1	4.1	15.1	10.8	2.9
	Permutation importance	1.4	51.6	18.0	2.0	4.6	14.8	6.9	0.6
RCP 2.6–70/ 83.3 ± 2.5	Percent contribution	6.8	52.5	11.4	5.2	6.9	5.5	9.6	2.1
	Permutation importance	4.8	35.2	30.8	6.9	8.9	6.6	6.1	0.7
RCP 4.5–70/ 80.8 ± 2.9	Percent contribution	10.3	28.9	4.3	32.2	5.7	3.8	11.9	2.8
	Permutation importance	2.4	39.5	9.7	26.0	4.6	12.5	4.8	0.6
RCP 6.0–70/ 79.3 ± 3.6	Percent contribution	9.8	31.9	31.9	34.9	2.5	1.3	12.6	3.3
	Permutation importance	5.4	57.4	57.4	3.9	4.0	4.9	9.4	1.0
RCP 8.5–70/ 81.0 ± 3.7	Percent contribution	11.2	51.2	4.1	14.3	1.8	3.3	11.6	2.5
	Permutation importance	2.0	36.3	9.4	31.8	1.1	14.6	4.4	0.4

). The average replicate runs were 0.836 ± 0.008, which revealed the best fitting ability of the model with the bioclimatic variables used for the test data. In addition, the habitat suitability range generated using a MaxEnt approach showed a moderate to very good correspondence for a predictive map of G. robusta, as revealed by the statistical measures, such as K (0.216), NMI (0.240), TSS (0.686), and others that are shown in

Table 3. Confusion matrix-derived measures of classification accuracy.

Measures	Calculated Value
Overall accuracy	0.822
Error rate	0.177
Prevalence (P)	0.037
Overall diagnostic power	0.963
Correct classification rate	0.822
Sensitivity (Sn)	0.865
Specificity (Sp)	0.821
False positive rate	0.217
False negative rate	0.135
Positive predictive power (PPP or TPR)	0.157
Negative predictive power (NPP or TNR)	0.994
Misclassification rate	0.178
Odds ratio	29.301
Kappa (K)	0.216
Normalized mutual information (NMI) n (s)	0.240
True skill statistic (TSS)	0.686

. The model revealed that the minimum temperature of coldest month (Bio 6), slope (Slo), mean temperature of driest quarter (Bio 9), and precipitation of driest quarter (Bio 17), contribute significantly and all the variables have associated permutational importance.

Geospatial test (used for validation) data were then overlaid on MaxEnt, which fell in the actual area of species occurrence, revealing a high level of accuracy in the prediction and estimation of distribution. The current habitat suitability map indicates that the maximum percentage of geo-coordinates that occur in the actual area of validation of the model mostly fell in the districts of Champawat and Tehri (100%), followed by Haridwar (80.0%) and Uttarkashi (78.04%) (Figure 2).

The MaxEnt modeling revealed the current habitat suitability area of ~1641 km² for G. robusta, which represents 3.07% of the total geographical area and 5.62% of the non-forest areas in the state of Uttarakhand (

Table 4. Estimated area under G. robusta in different districts of Uttarakhand Himalayas through MaxEnt modeling.

Sl. No.	Districts	Geographical Area (km2)	Non-Forest Cover (km2)	Estimated Area (km2)	Estimated Area (%) in Respect to Total Geographical Area	Estimated Area (%) in Respect to Non-Forest Cover
1.	Almora	3139	1426	317.4	5.87	1.09
2.	Bageshwar	2246	980	56.26	4.20	0.19
3.	Chamoli	8030	5321	113.97	15.01	0.39
4.	Champawat	1766	542	21.78	3.30	0.07
5.	Dehradun	3088	1483	257.56	5.77	0.88
6.	Haridwar	2360	1772	84.53	4.41	0.29
7.	Nainital	4251	1203	53.12	7.95	0.18
8.	Pauri	5230	1935	320.86	9.78	1.10
9.	Pithoragarh	7090	5012	52.38	13.26	0.18
10.	Rudraprayag	1984	843	87.51	3.71	0.30
11.	Tehri	3642	1577	50.55	6.81	0.17
12.	Udham Singh Nagar	2641	2106	61.83	4.94	0.21
13.	Uttarkashi	8016	4988	163.69	14.99	0.56
Total		53,483	29,188	1641.44	3.07	5.62

Source: Geographical area and non-forest cover, forest survey of India report (FSI 2019).

). The highest distribution area (~321 km²) was recorded for Pauri district, followed by Almora (~317 km²), Dehradun (~258 km²), and Uttarkashi (~164 km²), while the lowest distribution area was estimated for Champawat (~22 km²) and Bageshwar (~56 km²).

Furthermore, in order to investigate the possible relationships between species habitat suitability and climatic conditions, the current map was overlaid on the KGCC map. The KGCC map for the current distribution of G. robusta classified the Uttarakhand Himalayas into 5 subtypes (Figure 4). According to the generated output, the species occurs in four climate regions. Remarkably, the maximum occurrence of species was found in the humid subtropical climate (Cwa; C = warm temperate, w = winter dry, and a = hot summer) of the lower stretches of the Shivalik range, followed by the subtropical highland oceanic climate (Cwb; C = warm temperate, w = winter dry, and b = warm summer) of the middle Himalayas; however, a very low occurrence was observed in the monsoon-influenced warm, humid summer continental climate (Dwb; D = Snow, w = winter dry, and b = warm summer) and subarctic climate (Dwc; D = Snow, w = winter dry, and c = cool summer). This type of climatic condition prevails in the high-altitude ranges of the middle Himalayas of the districts Chamoli, Pithoragarh, and Rudraprayag.

3.3. Predicted Future (2050s and 2070s) Climatic Habitat Suitability Range

The future habitat suitability range of G. robusta is depicted in Supplementary Figures S2a–h and S3a–h. The environmental variable minimum temperature of the coldest month (Bio 6) contributes an overall maximum (52.8%; 51.6) during RCP 8.5–50, whereas the minimal (1.5%; 0.9) was attained by the precipitation of the driest month (Bio 14) during RCP 4.5–50. Compared with the area of the most optimal habitat, the predictions for 2050 and 2070 using the RCP 2.6, 4.5, 6.0, and 8.5 climatic models show almost or slightly more habitat suitability than the current prediction. The aerial enhancements of ~85 km² and ~226 km² were shown by RCP 2.6–2070 in 500 m buffer and RCP 6.0–2070 in 1000 m buffer, respectively, when compared with the current scenario.

Notably, under the future scenarios, in the 500 m buffer zone, the RCP 8.5 shows a maximum potential occurrence of ~109 km² and ~97 km² for the years 2050 and 2070, respectively, with predicted species distribution in the non-forest agricultural areas and wastelands in the hills (Figure 5a–d). The non-forest regions of Tons (Upper Yamuna, Barkot, Uttarkashi), Kalsi (Mussoorie, Dehradun), the Kedarnath Wildlife Sanctuary, and the Badrinath Forest Division are the potential areas for habitat suitability. For a 1000 m buffer, RCP 6.0 shows the maximum potential occurrence area (~309 km²) for 2050 and RCP 8.5 (~247 km²) for 2070. The species occurrence in non-forest, pasture, and other adjacent land areas is similar to the 500 m buffer zone with the inclusion of the Almora Forest Division during 2050.

3.4. Geological Correlation to Predict Current Distribution and Future Habitat Suitability Range for G. robusta

Remarkably, the map (Figure 6) shows that the Uttarakhand region is marked by a prominent NW–SE trend of geological variation, whereas the distribution of G. robusta rather follows a preferential NE–SW trend. The trend followed by the species coincides with the main orientation of the valleys carved by rivers that drain the more elevated areas. Importantly, the occurrence of G. robusta significantly coincides with soils formed over rock types with high iron content. Notably, the five main locations of the habitat suitability range could be identified. Location 1 consists of rocks from the Jaunsar Group (Mandhali, Chandpur, and Nagthat formations), which is composed of quartzite, slate, phyllite, dolomite, and meta-basic volcanics (high iron content). These are also influenced by rocks from the Central Crystallines, located upstream, which are composed of quartz–mica–chlorite–hornblende schist (high iron content) and gneiss. Location 2 comprises the Dum gravels and Kapkot alluvium, in which gravel, pebbles, and boulders are included in an oxidized matrix (the iron oxide-rich matrix makes iron promptly available). Location 3 consists of rocks from the Garhwal Group, limestone, slate, and quartzite associated with metabatic volcanics (purple color in the map, with very high iron content). Location 4 presents a similar context to location 3; it is in the Garhwal Group and has high iron meta volcanics. Lastly, location 5 is the weakest habitat suitability range, as it comprises rocks from the Almora Group and the Gorakhnath Formation, which includes garnet–mica–chlorite schist (high iron content, but not as high as the rocks at the other sites), quartzite, and phyllite.

4. Discussion

4.1. Geographical Distribution, Stand Structure, and Phenotypic Plasticity

The ecological niche of a species exhibits different environmental conditions within the distribution area; thus, it is not possible that only a single phenotype reveals the highest fitness levels (survival and adaptation) in all the conditions [67]. Species with wide geographical range have the potential to exhibit large intraspecific variations in morphology, which might be a good model to study the local and regional adaptations [68]. In India, the genus Eucalyptus, Populus, and Casuarina have been studied in greater detail with respect to their morphological traits and molecular variations, besides analyzing their adaptability to a particular climate, soil type, and region. No attempt has been conducted on the basic distribution and morphological attributes in G. robusta to reveal the species’ phenotypic plasticity.

Therefore, the present work on this species was carried out in Uttarakhand Himalayas, which consist of five climatic zones, i.e., tropical, montane subtropical, montane temperate, sub-alpine, and alpine, with an altitude range from 150 m to above 7800 m msl. Different morphological parameters, such as TH, GBH, CBH, CW, BA, DI, and IPI, were differentiated into four classes by assuming that the size class (DBH) is directly related to the age class distribution [69]. Under this assumption, the lower limit of the class interval was defined for juvenile individuals followed by young, adult, and mature trees. This approach might allow for more accurate classification, in addition to providing the health status of the stand, thus revealing the morphological trait plasticity of the sampled trees used in the study. The height and diameter class revealed that the young and adult trees were found at the maximum number in the mid-altitudinal range, i.e., 601–>1500 m msl (the maximum number of individuals falls within TH ranged between 9.1 and 27 m and DBH class from 0.24 to 0.44 m) with both parameters, i.e., TH and DBH, falling in the middle range class, showing the demographic stability of the species in the Himalayan terrain. Corroborating our findings, changes in the stand structure for TH and DBH were obtained for the vegetation present in the temperate forests of Mount Norikura, Japan [70].

Similarly, the effect of the DBH class distribution was observed in different altitudinal gradients in the western Himalayas [71]. Based on the number of individuals in different altitudinal classes, the study on the altitudinal behavior of the tree species (365 species) in the Subansiri district, the eastern Himalayas revealed that the highest frequency of species were restricted to 400–600 m and 800–1000 m msl [72]. In addition, the parameters, i.e., CBH, CW, CA, IPI, and DI, were also studied for G. robusta; however, limited studies have been carried out for variation in CBH, CW, and CA with respect to the altitudinal gradient, which also reveals that the maximum number of individuals found in the mid-class range and parameters falls mostly in the mid-altitudinal range (601–1500 m). Importantly, research has been conducted for DI and IPI in relation to the altitudinal gradient, where some species are resistant to a particular pathogen or insect at a lower altitude, while others are resistant at a higher altitude. For instance, infestation by Antestiopsis thunbergii in Coffea arabica was assessed in 24 coffee farms located at an elevation gradient of 1000–1700 m, where the bug prefers the host species at higher elevations [73].

4.2. Habitat Suitability Range through MaxEnt Modeling

Predicting the extent of the spread of an exotic plant in its introduced range is difficult due to a lack of detailed knowledge of the species’ biology and ecology, and the key environmental predictors in its new ranges; thus, an assessment of the habitat’s suitability is a valuable modeling approach used to predict the appropriate conditions for the introduction of any exotic tree species [74]. As G. robusta is mainly a plantation crop, prediction modeling and habitat suitability mapping for an appropriate cultivational/planted region through modeling-based approaches are now becoming an effective and applicable technique to assess the concurrence/relationship between environmental factors and the species [75].

The biomes are mainly controlled by a climate that strongly determines plant and animal distribution on the planet [76]; thus, climatic classifications provide a convenient tool for the validation of models for species distribution with a simulated analysis of future climate changes [77]. Notably, out of the total area (53,483 km²) of Uttarakhand, an area of ~1641 km² was estimated with a significant AUC value (0.836 ± 0.8) determined by the performance and parameters, such as K (0.216), NMI (0.240), and TSS (0.686), which revealed the accuracy of the MaxEnt model. In the present study, all the models had AUC values higher than 0.8, except RCP 6.0–70. Similar AUC values were reported for Vincetoxicum arnottianum [78], Hyptis suaveolens [79], Scutellaria baicalensis [80], and Oxytenanthera abyssinica [81]. Notably, the K value is the most widely used measure and shows fair agreement for the performance of models generating presence–absence predictions in G. robusta, whose statistical value ranges from −1 to +1, where +1 indicates perfect agreement and values of zero or less indicate a performance no better than random [82]. Despite its wide use, several studies have criticized it for being inherently dependent on prevalence, and researchers have argued that this dependency introduces statistical artefacts to estimates of predictive accuracy. Therefore, the NMI and TSS are better measures of the performance of presence–absence distribution models, which corrects for this dependency on prevalence while keeping all the advantages of K [66]. However, in a critique to Allouche, Somodi et al. (2017) [83] discussed that TSS is prevalence-independent, but only under strict assumptions and large sample sizes. A TSS value closer to one is considered excellent.

Determining the method by which the correct number of environmental factors is chosen is important for using MaxEnt to predict the distributions of exotic species. According to the contribution rate of environmental factors and the permutational importance in the jackknife test, this study found that the minimum temperature of the coldest month (Bio 6), the slope (Slo), the mean temperature of the driest quarter (Bio 9), and the precipitation of the driest quarter (Bio 17) predict the occurrence G. robusta in the Uttarakhand. In the case of Xanthium italicum, factors such as the annual mean temperature, the monthly mean diurnal temperature range, the standard deviation of the seasonal temperature change, and the annual average precipitation contributed most to the prediction of the exotic range of distribution in China [84]. Similarly, endemic Pomatosace filicula under climate change scenarios revealed that the altitude, annual precipitation (Bio 12), and annual temperature range (Bio 7) contributed most in determining the range of species in the Qinghai–Tibet Plateau [85]. Furthermore, the variables of temperature seasonality and isothermality contributed significantly in predicting the current and future distribution scenarios of Acacia reficiens and Opuntia spp. in Laikipia–Samburu, Kenya [86]. All these studies suggest that temperature variables contribute remarkably to the prediction of a species’ exotic range.

Though the G. robusta has a wider distribution range and, in our case, lies between 165 and 1957 m msl, the maximum probability of presence was limited to a range between 700 and 1600 m msl. Moreover, the species shows better adaptability across the altitudinal gradient (<601–>1500 m msl), which represents climatic conditions ranging from hot summers to cold winters. The resulting map effectively shows the current distribution in humid subtropical climate areas, particularly in the lower stretches and foothills of the Shivalik, followed by the subtropical highland oceanic climate areas of the middle Himalayas (Figure 4). In these elevation ranges, the precipitation varies from 200 to 250 cm per annum and the lower average temperature ranges between 19 and 21 °C. The lower montane zone has an average annual temperature of 14–18 °C with precipitation of 250–300 cm per annum [87]. These ideal climatic conditions were also revealed to match the widespread distribution of G. robusta in Rwanda [88]. The same classification system has been used to determine the global potential of exotics in different regions of the world, i.e., Jatropha curcas’ distribution, a plant that can be used to help replace fossil fuels [89]; and the system has also been applied when seeking potential distributions of native and alien Australian Acacias, i.e., A. cyclops and A. pycnantha in Australia and South Africa, respectively [90].

4.3. Future Habitat Suitability Range

Since the species is well adapted, as evidenced by its vast distribution and broad geographical range in the study area (Figure 2), it might be possible that the species shows the sign of potential suitability and adaptation in an extensive range of Uttarakhand and other regions of the northwestern Himalayas. Exotic tree plantations offer the potential to respond to increasing pressure for deforestation, besides delivering social, economic, and environmental advantages [91]. However, to understand the possible extent of the introduced species, it is necessary to ensure their broad range of establishment that may restore the wedged habitat. For such a study, habitat suitability modeling provides a tool for researchers and land managers to understand the possible ranges of introduction and make decisions about the future prospects of the species’ introduction [40]. The exotic species could be able to contribute to increased economic returns in addition to being capable of having a strong impact on the ecosystem processes, particularly when it differs from the native flora in many key morpho-functional features [92]. However, the adaptive behavior of a non-native species in a new ecological niche [93] realizes the lagging of natural pests or predators. The prevention of problems, such as species invasion, aggressive behavior, threat to native species, etc., can be achieved through the implementation of SDM, which are now increasingly used to predict the habitat suitability ranges of exotics [94]. Only a few model-based studies have been conducted for exotic tree species [95,96,97,98], but none of these were conducted in India.

Additionally, to overcome these challenges, the distinction between the predictions of current and future distributions is critical [94]. In the prediction modeling, the key environmental variables of habitat suitability included the climate (mostly represented by temperature seasonality; annual standard deviation in temperature), altitude, distance to the nearest urban and built-up area, etc. [99]. In order to map the suitability pattern, it is necessary to understand the factors controlling the species distribution [100]. This will help in managing exotics with less negative impact on the ecosystem. These facts suggest that the effective use of models in a risk assessment framework will greatly increase the utility of SDMs for managing exotic species in a judicious manner. Prominently, the predictive accuracy of SDMs was slight, but improved significantly when the models were calibrated with the habitat suitability data only, not considering the native species [100] (Barbet-Massin et al., 2018). The MaxEnt binary output provides evidence that a model is robust when it is fit with a limited dataset. It is an appropriate tool to use when land management objectives are supported by limited resources; thus, it requires a conservative but highly accurate estimate of habitat suitability for exotic species in the landscape [101,102].

Accurate monitoring of the distribution of any exotic species is of the utmost importance and is urgent and essential for developing cost-effective control strategies before the large-scale introduction of a species. Earlier research has also suggested the effective use of SDM and RS&GIS technology for depicting the potential habitat suitability pattern of different species around the globe. Many studies have been carried out for predicting the habitat suitability of natural forest tree species, such as Juniperus sp. in the southern Zagros Mountains of Iran [103], Perilla frutescens in Uttarakhand, India [104], and Carthamus tinctorius in China [75]; as well as for exotic tree species, such as Ligustrum lucidum in Córdoba, Argentina [105] and Robinia pseudoacacia in northeast Slovenia [106]; and the mapping of urban tree species, such as Acer platanoides, Thuja plicata, Prunus serrulata, Pseudotsuga menziesii, Gleditsia triacanthos, and A. rubrum in Columbia [107], and Spartina alterniflora in the central coast of Jiangsu in east China [107].

4.4. Geological Correlation to Current and Future Habitat Suitability Prediction

Finally, in addition to climate control, bedrock geology is also an important and significant factor that determines the soil composition, thus exerting control over the distribution and potential adaptability of any plant species [108]. As presented in Section 3.4, the results of the geological analysis reveal an intriguing geological correlation of the habitat suitability pattern of G. robusta with the occurrence of rocks presenting an elevated iron content. During the field surveys, fewer young seedlings were found, which might be due to the fact that the species mainly occurs adjacent to the major roads that are usually covered with concrete material. In some places, particularly in landslide-prone areas and soil dumping zones, some seedlings have been able to thrive. These observations demonstrate the species’ great potential to thrive in varied soil conditions; hence, the geological pattern of the bedrock in the region presents very low control over the distribution of G. robusta. The statement is also supported by the ability of G. robusta to grow in some of the poorest soils in the world, especially where phosphorus is limited [109]. Despite being able to grow in any soil type, G. robusta favors growing in fertile, deep, and light soils; however, the species shows poor growth in shallow, stony, and sandy soils, where the soil fertility is low [88].

Moreover, the species does not possess regeneration problems [110]. It has been known since the last century that some tree species of the Grevillea genus present intolerance to high phosphate availability, known as phosphate toxicity [111,112]. Increased calcium levels (in soils formed near limestone and dolomite rocks) can exacerbate toxicity, but sensitivity to high phosphate levels can be reduced by improving the iron supply, even when exacerbated by calcium [112,113,114]. According to these studies, an abundant iron supply can help the susceptible plants overcome some of the physiological problems associated with toxicity. This seems to be the case in our study area, as indicated by the potential adaptation prediction results. Figure 6 shows that the potential is more likely to occur in areas where the soil provides the means (iron) for the G. robusta to overcome phosphate toxicity, improving its growth, reproduction, and consequent success in adaptation. Calcium-rich rocks interlayered in most of the sedimentary sequences in the study area could be the cause of toxicity issues in this species. Another possibility is that G. robusta can use iron to acquire phosphate when the phosphate availability is extremely low; this is achieve by the mobilization of iron phosphates around its roots, a process that was demonstrated experimentally [115].

5. Conclusions

This study provides novel information on the size-class distribution of G. robusta plants, which were found to vary over time depending on the biological and environmental conditions. Importantly, the frequency distribution of the diameter class ranged from 0.24 to 0.44 m, indicating that young and adult trees were more recurrent. This shows the demographic stability of the species in the northwestern Himalayas. This study on G. robusta also provides novel information on the current distribution pattern, a genetic resource assessment, future prediction scenario, and potential habitat suitability range in the Himalayas through MaxEnt modeling approach. Under this scenario, the KGCC suggests humid subtropical climate (Cwa) conditions prevail for the maximum occurrence of the species; however, the natural forest areas of Tons (Upper Yamuna Barkot, Uttarkashi), Kalsi (Mussoorie, Dehradun), the Kedarnath Wildlife Sanctuary, Almora, and the Badrinath Forest Division show the potential habitat of G. robusta in different RCPs for the years 2050 and 2070. This distribution pattern of the species was overlaid and bestows low correlation to the geological pattern of the bedrock in the region, revealing the ability of G. robusta to grow in any soil type. The habitat suitability ranges are mostly confined to the regions with the occurrence of rocks (and soils) with an elevated iron content. Iron is possibly used by G. robusta as a condition to overcome sensitivity to toxicity, a common issue in the Grevillea genus that might be amplified by the presence of calcium-rich soils that are also present in the region. Iron can also be used by this species for the mobilization of iron phosphates in extremely low P soils. Irrefutably, the present and future habitat suitability maps, along with influential bioclimatic variables could be of paramount importance to the forest department and land managers, policy makers, and research and development organizations in India and around the world, who may use the results from this study for planning long-term sustainable forest management policies for exotic tree species.