Determining bioclimatic space of Himalayan alder for agroforestry systems in Nepal.

Himalayan alder species are proven to be very useful in traditional as well as contemporary agroforestry practice. These nitrogen-fixing trees are also useful in the land restoration. Therefore, understanding the distribution of Himalayan alder and the potential zone for plantation is meaningful in the agroforestry sector. Suitable climatic zones of Alnus spp. were modelled in MaxEnt software using a subset of least correlated bioclimatic variables for current conditions (1950-2000), topographic variables (DEM derived) and Landuse Landcover (LULC) data. We generated several models and selected the best model against random models using ANOVA and t-test. The environmental variables that best explained the current distribution of the species were identified and used to project into the future. For future projections, ensemble scenarios of climate change projection derived from the results of 19 Earth System Models (ESM) were used. Our model revealed that the most favorable conditions for Alnus nepalensis are in central Nepal in the moist north-west facing slope, whereas for Alnus nitida they are in western Nepal. The major climatic factor that contributes to Alnus species distribution in Nepal appears to be precipitation during the warmest quarter for A. nepalensis and precipitation during the driest quarter for A. nitida. Future projections revealed changes in the probability distribution of these species, as well as where they need conservation and where they can be planted. Also, our model predicts that the distribution of Alnus spp. in hilly regions will remain unchanged, and therefore may represent sites that can be used to revitalize traditional agroforestry systems and extract source material for land restoration.

Introduction

The importance of carbon sequestration increases when there is a sufficient number of trees to absorb the emitted carbon (McCarl and Callaway, 1993). Absorption of carbon is only possible when nitrogen binds phytomass with assimilated carbon during photosynthesis. Fixation of atmospheric nitrogen is essential so that it is readily available to trees (Galloway et al., 2004). The genus Alnus, belong to the birch family host species from the symbiotic endophytic genus Frankia (Khan et al., 2007). The symbiont Frankia uses carbohydrates from alder trees to convert atmospheric N2 into reactive nitrogen, a nitrate form directly available to plants growing in the soil (Myrold and Huss-Dannel, 1994). Alnus species host nitrogen-fixing bacteria; therefore, they are ecologically important in nitrogen fixation. Alnus are pioneer species that establish in eroded and exposed soils, which in the course of time improve land health. Alders are vigorous and fast-growing, even in acidic soil and damaged sites such as burned areas and mining sites, further adding to their importance as the species of choice in forest restoration programs.

Alnus spp., specifically A. nepalensis, have commonly been used in traditional agroforestry systems as shade, fodder, fuelwood and timber. Traditional agroforestry systems are unique in hilly regions of Nepal and adjoining Himalayan countries (Sharma et al., 2007). Previous research has revealed that this species is an excellent choice for soil restoration and sloping land management. Alnus plantation with cardamom (Sharma et al., 2002), tea (Mortimer et al., 2015), and mandarin oranges (Chhetri and Gauchan, 2008) represent a few examples of the importance of Alnus in the agroforestry system. Such mixed plantation is beneficial to commercial crops.

Agroforestry systems are examples of the ecological aid and economic viability for people compared to rainfed agriculture (Sharma et al., 2007). Therefore, it is crucial to discover suitable areas for important tree species such as Alnus for agroforestry planning (Ranjitkar et al., 2016). Understanding and identifying potential growing areas will provide insight for the intercropping of commercial plant species that are essential for the livelihood of smallholder farmers. Alder species are planted throughout the hilly region and distributed throughout the Himalayas. It is essential to accurately discern how climatic factors determine the distribution probability of different Alnus species to recommend right species at the appropriate bioclimatic zone. Also, bioclimatic zone identification can identify suitable spaces where these species can be grown alternative to each other. In this study, we aim to describe the habitat suitability of the alder using species distribution modelling (SDM) methods. Specifically, we determine the best model for distribution probability and habitat suitability of focal species. We use that model to identify the climatic variables critical for occurrence of Alnus species. Finally, we determine the areas for plantations of two Alnus spp. in Nepal for appropriate agroforestry practice.

Materials and methods

Focal species

Two species of Alnus are reported in Nepal: Alnus nepalensis D. Don and Alnus nitida (Spach) Endl. A. nepalensis is the eastern Himalayan component distributed from northwest India through the Himalayas to Guangxi, Guizhou, Sichuan in China in the east, to Thailand, Vietnam, and Laos in the south (Shaw et al., 2014a). A. nitida is a native to the western Himalayas distributed from Afghanistan and Pakistan to northwestern India and western Nepal (Shaw et al., 2014b). In Nepal, A. nitida is somewhat rare and still requires careful investigation in the western part of the country. The species is known to occur in pure forest stands in the upper Karnali river valley in western Nepal (Fig. 1). It grows well in moist sandy soils, whereas A. nepalensis grows throughout hilly regions, in a wide range of soil including freshly exposed landslide, as well as degraded and disturbed sites. Encroachment, forest degradation, over-exploitation for fuel-wood and timber collection are known threats to A. nepalensis trees growing in the Nepal. Biological threats to A. nepalensis include the stem borers Batocera spp., an aphid Eutfichosiphym alnifoliae pest of economic importance, and a parasitic plant infestation (O'Neill and Rana, 2016). Biological threats have yet to be identified for A. nitida; however, forest degradation and habitat fragmentation are known to threaten the species.

Fig. 1

Himalayan alder records from different parts of Nepal. Vertical lines represent tentative division of eastern (east to 86.3°), central and western (west to 83°) Nepal based on phytogeography (Stearn, 1960).

Species data

We used herbarium specimen data (from 1968 to 2012) and field surveys (Appendix A4) to compile the location of occurrence for the focal species. Specifically, herbarium specimen data were gathered directly at the National Herbarium and Plant Laboratory (KATH) and Tribhuvan University Central Herbarium (TUCH) (Appendix A1, Appendix A5). In addition, specimen data were gathered online from University of Tokyo (TI) (Appendix A2) and the Royal Botanical Garden at Edinburgh (RBGE) (Appendix A3, Appendix A6). Present geographical locations of the species were confirmed through “ground truthing” in several places, where rapid field surveys, communication with community forest members, forest officers and residents were useful. All the geo-referenced recorded points were used to check spatial autocorrelation of the species distribution using binary presence/absence data. We selected random background points within an area of 50 km radius from the point of occurrence for each focal species in order to examine Moran's I for spatial autocorrelation in Arc GIS 10.2.1. Presence points were spread to 10-min spatial grid to obtain least spatial autocorrelation (Ranjitkar et al., 2014a), which reduces A. nepalensis occurrence to 67 and A. nitida to 12 (Appendix B1). We use these treated points for further analysis in the selection of predictive variables and the calibration of the model.

Environmental and bioclimatic data

To model the potential distribution of Himalayan alders, we first collected 26 predictive variables (21 bioclimatic, four topographic and LULC). Out of 26 variables, we selected significant variables for modelling, considering their contribution and importance in the model calibration, collinearity among variables, Jackknife analysis and biological significance for the focal species. Among the two highly correlated variables (Pearson correlation coefficient, r > 0.80) (Appendix C1), one was chosen for its biological relevance to the species and ease of interpretation for inclusion in the model (Kumar and Stohlgren, 2009, Rana et al., 2017). Finally, a subset of the least correlated variables was selected (Table 1) as predictor variables.

Table 1

Variables used as inputs to generate distribution models for alder species in Nepal.

Category	Name of variables	Source
Years	Current as baseline and 2050 for future
Scenarios	Ensemble of downscaled CIMP5 ESM Models under 4 different RCPs	Taylor et al., 2012
Area	Nepal (1,47,181 km2)
Spatial resolution	1 km (∼30 arc-sec)
Climate	Bio3 = Isothermality	www.worldclim.org/bioclim
	Bio4 = Temperature Seasonality	Hijmans et al., 2005
	Bio6 = Min Temperature of Coldest Month
	Bio9 = Mean Temperature of Driest Quarter
	Bio15 = Precipitation Seasonality
	Bio17 = Precipitation of Driest Quarter
	Bio18 = Precipitation of Warmest Quarter
Topographic	Slope	www.eros.usgs.gov
	Aspects (Northness)	Calculated from DEM
	Aspects (Eastness)
	t_rgh = Terrain roughness
LULC	LULC = Landuse Landcover	www.diva-gis.org/gdata (GLC 2000)

The bioclimatic dataset available from the Worldclim (www.worldclim.org/bioclim; Hijmans et al., 2005) was used for current (1950–2000) predictions. At the same time, future (2050) bioclimatic variables were derived from the results of 19 Earth System Models (ESM) provided by the Coupled Model Inter-comparison Project-Phase5 (CIMP5; Taylor et al., 2012). Slope, aspects, and terrain roughness derived from Aster Digital Elevation Model (DEM) (www.eros.usgs.gov) and Landuse-landcover (LULC) (Roy et al., 2003) were also included in the model (Appendix C2). All the predictor layers are at a spatial resolution of 30 arc-seconds.

Modeling approach

We used MaxEnt 3.3.3k (Maximum Entropy) for modelling purposes. MaxEnt uses only presence data in combination with environmental data for the whole study area to derive a model and predict suitable conditions or ecological niches (Phillips et al., 2004, Phillips et al., 2006, Phillips and Dudik, 2008). Altogether 50 replicates of the model with the least-correlated set of variables were run, setting MaxEnt to select 75% of occurrence localities randomly for calibrating, and the remaining 25% for evaluating at each run. The bootstrap method was used with maximum iteration = 5000, convergence threshold = 0.0001, maximum background points = 10,000, replication 50 and default features.

The current suitable niche was forecasted using the significant models among the four models using least correlated variables. The four models used bioclimatic variables, bioclimatic with topographic, bioclimatic with LULC and bioclimatic with topographic and LULC.

Future projection

The potential future distribution probability and habitat suitability of focal species was forecasted using an ensemble scenario of climate change projection derived from the results of 19 Earth System Models (ESM) provided by the Coupled Model Inter-comparison Project-Phase5 (CIMP5; Taylor et al., 2012). There are four representative concentration pathways (RCP; Vuuren et al., 2011) within each of the 19 ESM, ranging from RCP 2.6 (aggressive mitigation/lowest emissions) to RCP 8.5 (higher emission scenario). All models within each RCP were ensemble using the class with the majority of occurrence within any particular grid cell as the class for location (Zomer et al., 2014, Ranjitkar et al., 2016). For the modelling process, LULC layer was considered as a stable, as the projected for the future is not available (Ranjitkar et al., 2016). It is necessary to measure the similarity between the new climatic conditions and those in the calibration to make a projection for places and time not sampled in calibration data (Elith et al., 2010). The new predicted environment is an extrapolation. The extent of extrapolation in the future projection was examined using multivariate environmental similarity surface (MESS) (Elith and Leathwick, 2009, Elith et al., 2010).

Spatial analysis

We determined spatial changes in the habitat suitability of focal species by transforming the probability of occurrences from models into presence-absence maps. The resulting raster layers were reclassified as ‘1’ for the presence and ‘0’ for absence using binary maps based on a threshold. The logistic output in Maxent modelling with ten percentile presence was a threshold value in this work. The threshold represents the cut-off point above which climate is favourable for the focal species. We reclassified the MaxEnt output rasters into three classes of habitat suitability: low suitability (10–50% probability of occurrence), medium suitability (50–75% probability of occurrence) and high suitability (>75% probability of occurrence). We differentiated all the future predictions into stable, retracted and expanded areas, comparing them to current predictions. All the spatial analysis was performed in ArcMap 10.2.1 (ESRI) and the extension Spatial Analysis (ESRI).

Statistical test of significance of model

The evaluation of the accuracies of the prediction models generated through MaxEnt was carried out selecting the threshold-independent Receiver Operating Characteristics (ROC) Area Under Curve (AUC) method (Liu et al., 2005, Phillips et al., 2008). The AUC values from the 50 subsets produced by each model were statistically tested to determine if they were significantly better than random AUC (0.5). The AUC averages were calculated and compared with different models.

We examined the mean AUCs of each model against the random model (AUC = 0.5), and other individual models generated. Based on the normality test, a one-tailed T-test was used to test for the significance of the average AUC values against a random model. Analysis of Variance (ANOVA) was conducted to check any significant difference between the means of the AUC. We selected the best model based on the pair-wise comparison with the random model and remaining models.

We formulated two hypotheses:

Testing the hypothesis that the models produced for each species are significantly better than a random model, i.e. H0: AUC (train) and AUC (test) = 0.5; H1: AUC (train) and AUC (test) > 0.5.

Testing the hypothesis to compare the means of the models, i.e., H0: AUC1 = AUC2; H1: AUC1 ≠ AUC2 (where, AUC1 = models with climatic variables and AUC2 = model that include topographic variables)

Results

The spatial autocorrelation performed using Moran's I for two species of Alnus in the form of binary prediction indicates that the points are randomly distributed. The Moran's I value for both species are less than 1 (−0.00521 for A. nitida; 0.022107 for A. nepalensis; P > 0.05) indicating negligible spatial autocorrelation. The Shapiro–Wilk test revealed a P-value greater than 0.05 except for the Test AUC1 (0.008) and Test AUC4 (1.281 e−6) of A. nitida (Appendix C3), indicating the training AUC for both the species distributed normally.

Threshold independent evaluation of the models

We found, the distribution probability and habitat suitability of A. nepalensis and A. nitida can be predicted using bioclimatic variables with topographic variables and LULC to achieve both test and training AUC that are significantly better than a random model. Thus, the null hypothesis H0 was rejected and the alternative hypothesis H1 accepted. The P-values calculated on the average of training and test AUC for all four models were significantly better than the random model (P < 0.00001, one sample T-test, 95% CI). The training AUC and training gain were highest for Model4 for both species (Table 2). Thus, Model4 (see Table 2, Table 3) was selected as the best model defining distribution probability and habitat suitability of Himalayan alders.

Table 2

Results of threshold independent evaluation and P-values of average AUC.

Species	Model	Training AUC	Test AUC	Training gain	Test gain	P-values of Average AUC
Alnus nepalensis	1	0.830	0.767	0.613	0.515	<2.2e−16
	2	0.868	0.780	0.780	0.543	<2.2e−16
	3	0.840	0.773	0.679	0.536	<2.2e−16
	4	0.889	0.772	0.959	0.469	<2.2e−16
Alnus nitida	1	0.879	0.852	0.707	0.847	<2.2e−16
	2	0.908	0.842	0.840	0.830	<2.2e−16
	3	0.869	0.827	0.668	0.754	<2.2e−16
	4	0.946	0.868	1.197	1.146	<2.2e−16

Table 3

Results of the pair-wise mean comparison of the four models developed per species (P-values shown).

Species		Pair-wise mean comparison of the models
		1,2	1,3	1,4	2,3	2,4	3,4
Alnus nepalensis	AUCtr	<2.2e−16	0.00972	<2.2e−16	1.34e−10	2.031e−07	<2.2e−16
	AUCts	0.1036	0.4834	0.5887	0.4298	0.3767	0.9033
Alnus nitida	AUCtr	0.0001582	0.2648	<2.2e−16	6.10e−06	7.638e−10	2.662e−16
	AUCts	0.5357	0.09713	0.4209	0.3558	0.2317	0.05571

Note: AUCtr = training AUC and AUCts = test AUC.

Sensitivity analysis

Jackknife analysis showed that the environmental variables with the highest gain, when used in isolation, are bio9 and bio6 for A. nepalensis and slope for A. nitida, which therefore retains the most useful information. For both species, the environmental variable that decreases the model predictability gain when omitted is the slope, which therefore appears to have the most information that is not present in the other variables (Appendix B2). The percent contribution (Fig. 2) estimated that bio18 and slope contributed highly in the model that determines suitable habitat for both Alnus species in the model. Similarly, LULC was an equally important limiting factor for both species in the model.

Fig. 2

Percent contribution for each variable calculated from model 4 produced for both the focal species (t_rgh: terrain roughness; LULC: Land Use Land Cover).

The response curves (Appendix B3) showed that the probability of A. nepalensis occurrence increases with the increase in bio18 (optimum between 500 and 1000 mm) (Fig. 3a) as well as bio3 with 47% temperature variability. Bio4 either larger or smaller than optimum value 4000 decreases the presence probability (Fig. 3c), whereas bio6, bio9, and bio15 have the optimum value 2.5 °C, 11 °C and 100 respective units showing the hump-shaped curve. LULC seems to have fluctuation in the probability of presence.

Fig. 3

Variable response curves for the top two predictors for Alnus nepalensis (a, c) and Alnus nitida (b, d). (Precipitation (Bio17 and Bio18) is expressed in mm and temperature seasonality (Bio4) in standard deviation*100).

In the case of A. nitida (Appendix B4), the habitat suitability increases with an increase in bio17 with rainfall of at least 80 mm (Fig. 3b) as well as terrain roughness with an optimum value about 1. For the remaining variables, like bio18 (Fig. 3d), bio3, bio4, bio6, bio9, bio15, eastness, and northness, habitat suitability decreases with an increase in their corresponding value. Slope around 23–27° was most suitable for this species. Similarly, habitat suitability increased in the needle-leaved evergreen forest, with closed to open tree cover.

Comparison of the means of the models

ANOVA results showed a significant difference in the means of AUC of all four models for A. nepalensis (F-value = 19.03, P < 0.05) and A. nitida (F-value = 12.3522, P < 0.05). Further, a pairwise comparison of the average values of the models (Table 3) showed the training AUC for all pairs of means of models produced for A. nepalensis were significantly different (P < 0.05) from each other. In the case of A. nitida, the pair of the average of model 1 and model 3 did not show any significance for training AUC. Similarly, test AUC for all pairs of the average of models were also insignificant.

Contrary to the training AUC, a comparison of the means of the test AUC did not show any significant difference for all model combinations for Alnus species. Therefore, based on the pairwise comparison of the means of the models, it can be concluded that all four models have the same predictive power regarding the test AUC for Alnus species (the test means were not significantly different from each other).

Gains produced from MaxEnt modelling have also been used in determining the best performing model (Yost et al., 2008). Gains describe how well the model fits the training and test dataset available. The average gains for all four models for both species of Alnus showed that model 4 achieved the highest average training gains. The AUC for A. nepalensis was 0.889 and for A. nitida 0.946 in model 4 (Table 2).

Current suitability

The model predicts clear climatic space for two species. East to 83.5° longitude was highly suitable for A. nepalensis, although habitat suitability extended toward the west as well, whereas A. nitida was confined west to 84° longitude. The climatic spaces for both species were along the river valleys. The suitable area for A. nepalensis was predicted to occupy 24% of the total area of Nepal in the baseline climatic scenario (Table 4). The central region of Nepal contains the largest and most continuous area for the species, whereas the western region has a dispersed distribution pattern. Furthermore, out of the total distribution probability and habitat suitable area, only 0.87% has high suitability and is mostly distributed in the central/eastern/upper hill regions; in contrast, 91.58% of the low suitability areas are distributed in the western/lower hill region (Fig. 4a, Table 4). Several isolated patches of varying but relatively small size occur mostly in the west and to some extent in the Siwalik range (Fig. 4a; c). Our model predicts the Terai and High Mountains as unsuitable areas for A. nepalensis.

Table 4

The projected areas of Alnus nepalensis and Alnus nitida under current and future (2050) climate change scenarios of CIMP5, derived from 19 ESM (four RCPs) categorizing low, medium and high suitability value.

Suitability	Low (10%–50%)		Medium (50%−75%)		High (>75%)
	km2	%	km2	%	km2	%
Alnus nepalensis
Current	32936	91.58	2713.7	7.55	313.41	0.87
RCP 2.6	22081	91.74	1720	7.15	268.96	1.12
RCP 4.5	21498	92.7	1520.3	6.56	172.52	0.74
RCP 6.0	20208	91.62	1632.6	7.4	215.47	0.98
RCP 8.5	19699	91.95	1533.9	7.16	189.85	0.89
Alnus nitida
Current	3748.8	88.4	423.4	9.98	68.56	1.62
RCP 2.6	9005.9	94.01	505.52	5.28	68.56	0.72
RCP 4.5	5619.5	88.29	614.76	9.66	130.33	2.05
RCP 6.0	3399.3	91.45	266.7	7.17	51.23	1.38
RCP 8.5	7624.2	90.88	672.77	8.02	92.67	1.1

Fig. 4

Predicted potential distribution of a) Alnus nepalensis and b) Alnus nitida under current bioclimatic conditions and c) combined habitat suitability of focal species.

Mountain of Western and High Mountain in central Nepal were climatically suitable for A. nitida (Fig. 4b; c). The total suitable area covers 0.03% of the total modelled area (Table 4). The western region was highly suitable for A. nitida, covering an area of 68.56 km2 (1.62%), whereas areas of low suitability were confined to 3748.82 km2 (88.40%) of the total suitable area (Table 4). The climatically suitable areas were predicted to be more within the extent of the Himalayas in the hilly districts in the Karnali and Dhaulagiri (Fig. 4b; c). The lower high mountains of central region are potentially a novel location for A. nitida, as there were no official reports or collection of the species from the area. The model also predicts some ecotones in the eastern hilly region as new places that lack the actual distribution of this species.

Future potentials

The model indicates that current habitat suitability for Alnus species will change greatly in the future (Fig. 5). The shifts in the distribution probability and habitat suitability vary among two species. A. nepalensis was predicted to move towards the northwest in the high mountains under all four RCPs (Fig. 5a–d). Under an ensemble scenario of climate change, A. nepalensis gained maximum suitable area in RCP 2.6 (about 24,070 km2) and declined to the lowest of about 21,423 km2 in RCP 8.5 (Table 4). The decline of habitat suitability of A. nepalensis was predicted to be confined mostly to central Nepal compared to eastern and western Nepal. The high and medium suitable areas were predicted to remain in the eastern and central regions, whereas a small fraction of medium and high suitability areas was projected towards the western region under RCP 2.6 and 6.0 (Fig. 5a–d). The habitat suitability of A. nitida was projected to likely increase (0.07% of the total area) by 2050 toward southeast under RCP 2.6 (Fig. 5e, Table 4). The projected suitability decreased under RCP 6.0 by 0.03% of the total suitability area (Fig. 5g, Table 4), whereas under other RCPs, the projected suitable area increased. The habitat suitability of A. nitida was predicted to occur more towards northwest districts of Nepal (Fig. 5e–h).

Fig. 5

Predicted future distribution of Alnus nepalensis (a–d) and Alnus nitida (e–h) in future climate scenarios of CIMP5 derived from 19 ESM under different RCPs.

According to the model, only 13.59% and 2.52% of the total area of Nepal under RCP 2.6 can be considered stable through 2050 for A. nepalensis and A. nitida respectively (Table 5). By 2050, our model (RCP 8.5) projected more reduction (13.77%), more expansion (3.83%) and less stability by 10.66% in the total area of Nepal compared to the current situation for A. nepalensis (Table 5, Fig. 6a–d). A. nitida was projected to show more expansion (3.97%) under RCP 2.6, reduction (1.17%) under RCP 6.0 and stability (2.52%) under RCP 2.6 (Table 5, Fig. 6e–h).

Table 5

The Retracted, Stable and Expanded areas for potential distribution of Alnus nepalensis and Alnus nitida by 2050 under four RCPs.

	Reduction		Expansion		Stable
	km2	%a	km2	%a	km2	%a
Alnus nepalensis
Current RCP 2.6	15964.88	10.85	4011.75	2.73	19996.97	13.59
Current RCP 4.5	17447.53	11.85	4593.36	3.12	18514.31	12.58
Current RCP 6.0	17621.56	11.97	3668.21	2.49	18338.78	12.46
Current RCP 8.5	20268.18	13.77	5630.76	3.83	15693.66	10.66
Alnus nitida
Current RCP 2.6	531.13	0.36	5837.19	3.97	3708.89	2.52
Current RCP 4.5	1233.28	0.84	3349.53	2.28	3006.74	2.04
Current RCP 6.0	1725.24	1.17	1194.11	0.81	2515.54	1.71
Current RCP 8.5	907.82	0.62	5037.1	3.42	3332.95	2.26

Percentage of reduction, expansion and stable area is calculated out of the total area of Nepal i.e., 147181 km2.

Fig. 6

Maps of the potentially stable, retracted and expanded areas with respect to the current distribution of focal species (a–d: Alnus nepalensis; e–h: Alnus nitida) under different RCPs (The future bioclimatic variables were derived from the results of 19 Earth System Models (ESM) provided by CIMP5).

MESS analysis (Fig. 7a–d) indicated some of the extrapolated regions along the western and central hilly region and in the Siwalik range will experience novel climates by 2050. The greater negative MESS values for the Terai region indicated points outside the training range, whereas positive MESS value in the hilly area indicated locations within the training range. The extrapolation is mostly along the river system.

Fig. 7

MESS analysis indicating novel environment across different RCPs: (a) RCP 2.6, (b) RCP 4.5, (c) RCP 6.0 and (d) RCP 8.5.

Discussion

In this study, we used a species distribution modelling (SDM) technique not just to understand the species environmental requirements but also to understand aspects of biogeography of Alnus species in Nepal. SDM is widely used in regions where there are least known occurrences of species (Elith et al., 2006, Guisan et al., 2006, Elith and Leathwick, 2009, Ranjitkar et al., 2014a, Zhang et al., 2012). The uncertainties and limitations in SDM techniques lead to under- or over-predicting species distribution (Barry and Elith, 2006, Ranjitkar et al., 2014a), and certainly create obstacles to generating accurate models. While gathering presence only data, spatial autocorrelation could have an impact on SDM (Ranjitkar et al., 2014a). The spatially correlated points lead to predicting the same habitat. Multiple observations within a confined geographic area could spatially correlate, as exemplified by A. nitida in Western Nepal. More importantly, when selecting background points within the region of interest for geographically restricted species, spatial autocorrelation should be prevented, which is properly addressed in the present work.

Physiological tolerance to climatic factors can influence distributions of plants (Woodward, 1992). Precipitation of the warmest quarter (for A. nepalensis) and precipitation of the driest quarter (for A. nitida) separates climatic spaces for Himalayan Alders. A. nepalensis found in the wetter part of Himalayas may require more rainfall than A. nitida found in the drier western Himalayas. Therefore, rainfall during the warmer and drier parts of the year are the important determining factors for distribution probability and habitat suitability of A. nepalensis and A. nitida respectively (Rana, 2014). This is the first study to investigate the relative importance of precipitation and temperature in order to determine Alnus species niches in Nepal. Areas that receive precipitation (during warmest quarter of the year) between 400 and 900 mm are good for occurrence of Himalayan alders. Sakalli (2013) predicted the global distribution of Alnus species was limited to areas of 500–1100 mm annual precipitation and temperatures between 8 and 15 °C.

The habitat suitability might not necessarily be an indication of species occurrence, as they are limited by edaphic and biogeographic factors (Ranjitkar et al., 2014b). The actual distribution is prone to limitation by seed dispersal mechanism, colonization capacity and LULC change (Sexton et al., 2009, Ranjitkar et al., 2014b). Extensive land use conversion in Nepal is one of the major reasons for the destruction of natural forests that impact the distribution of the least-explored and narrowly-distributed species like A. nitida. Fragmented habitat due to land conversion leads to a clustered distribution of this species, indicating the need to conserve the natural habitat. Mountainous topography and a river valley clearly divide the potential geographic region of Alnus species, possibly limiting seed dispersal and colonization. Climate change might limit or expand habitat suitability for the species.

Future projections indicate significant changes by 2050 in the habitat suitability of both Alnus species. The geographical origin of the genus Alnus and the role of climate in shaping their distribution is uncertain, but our model revealed expansions and reductions in suitability in response to climate. In addition to climate, high mountains (>3000 m asl) in the Northern aspect and increasing elevation prevent A. nepalensis seed dispersal and establishment, whereas an elevational shift of A. nitida might lead to the summit trap phenomena (Salick et al., 2009). The cold-adapted A. nitida, stressed by climate warming, migrate upwards until they reach the highest elevations, where they compete with other species. Sakalli (2017) predicted a northward shift of Alnus spp., and Ranjitkar et al. (2016) projected a contraction in the lower limit of A. nepalensis, both of which provide further support for our results.

The stable area predicted in our model for both species by 2050 can be designated areas for conservation of the species and as a source of germplasm, whereas new suitable areas can be designated zones for promoting an alder tree-based agroforestry system (Rana, 2014). The model may be essential for adaptive responses and preventive measures for the sustainable management of agroforestry tree species in the future. Agroforestry systems help farmers adapt to climate change by mitigating its consequences through carbon sequestration (Luedeling et al., 2014). Agroforestry systems are more complicated than monoculture situations. They consist of annual and perennial plants, which are often integrated with livestock (Luedeling et al., 2014). With respect to Nepal, these multipurpose agroforestry tree species are beneficial for large cardamom, which is an important perennial cash crop. Cardamom requires shade, therefore, intercropped under various types of trees is beneficial. A. nepalensis is the chief shade-providing trees used in large cardamom-based traditional agroforestry systems in the region (Sharma et al., 1998). However, the shifting trends in habitat suitability of A. nepalensis might cause depletion in the existing agroforestry system in Nepal. Insects defoliate most of the seedlings; and recently O'Neill and Rana (2016) reported parasitic (Loranthaceae group) infection of an A. nepalensis population in the hilly regions in central Nepal. Such infestation has not previously been documented in Nepal. Therefore, according to projections from our model, A. nepalensis populations in the stable region should be the focus of management. Also, A. nitida can act as a suitable substitute shade- and nitrogen-fixing tree for an intercropping Alder–Cardamom, Alder–Tea system, particularly in central Nepal where the climatic suitability for A. nitida is projected by 2050 in our model. It is the most suitable farm forestry tree in the high valley and mountains because of its frost hardiness, fast growth and nitrogen-fixing properties (Sheikh, 1993). A. nitida can be easily reproduced from seed and by vegetative means. It has no known insect or pest problems.

The Department of Forestry, the government of Nepal, have compared the growth rate of A. nitida with that of A. nepalensis in eastern Nepal (based on our personal communication). They reported 90% of seedlings of A. nitida survive attaining a mean height of 87 cm after two years. The growth performance of A. nitida closely matches or even exceeds that of A. nepalensis. According to our model, A. nepalensis will decline in the future, whereas A. nitida will gain suitability toward the north-east. Therefore, A. nitida could be a potential substitute for A. nepalensis at some high elevation places where parasite problems are reported. We recommend checking compatibility of this species with cardamom and tea, and establishing this rare species as a nitrogen-fixing shade tree for agroforestry systems in the future.

Conclusion

Our rigorous statistical model was successful in delineating the niche of focal species. Probably, A. nepalensis found in wetter parts of the Himalayas require more rainfall than A. nitida found in the drier western Himalayas. Therefore, rainfall during warmer parts and drier parts of the year are determining factors for habitat suitability for A. nepalensis and A. nitida respectively. Thus, wide variation in the species is due to precipitation, where A. nitida receives much less rainfall throughout the year (optimum 70 mm).

The uniqueness of this study is describing the suitable growing areas for Himalayan alders. Thus, the maps generated from the model should be considered as baseline information for the focal species. The most accurate predictive map for each species was used to identify potentially suitable niches as the core area. The results can be valuable inputs for prioritization of the areas for protecting the natural population and revitalizing the alder-based agroforestry system.

Collection date	Locality	District	Altitude (m)	Collection no.	Collector(s)
8/17/1997	Bhotebas–Gogane–Chichila–Kaptane	Sankhuwasawa		9755094	S. Noshiro et al.
2/12/2003	Palung–Daman, Nepal		1056		T. Watanabe et al.
10/29/1995	KhaniDanda–SheraKhalle–Diktel–MiyaKhola–DorpaChuriDanda	Khotang	1450	9555078	M. Mikage et al.
9/21/1996	RituKharka–Lulang–Lumsung–below Lumsung	Myagdi	2180	9686196	M. Mikage et al.
9/1/1996	Ratnechour–Beni	Baglung	830	9686008	M. Mikage et al.
1/19/1996	Sauraha–Padampur–Bhawanipur–HardaKhola	Chitawan		9611796	M. Mikage et al.
9/2/1995	Junbesi–Tagtor – a pass–Lanjura Pass–Goyom–Dakchu–Sete–Kensa	Solukhumbu		9596597	F. Miyamoto et al.
8/29/1971	near Nigale		1600	1050	M. Ohsawa & P.R. Shakya
10/23/1971	above Nhorbu Gaon		2100	997	M. Ohsawa & P.R. Shakya
5/5/1976	Royal Chitwan National Park	Chitawan	390		R.G. Troth
9/23/2007	Simigaun (2030 m)–Tamakoshi bridge (1430 m)–Lamabagar (1970 m)–Gungurkhola (1300 m)	Dolakha	2110	20740712	S. Noshiro et al.
9/13/2007	Dolakha (1750 m)–Nagdaha (890 m)–Malephu (940 m)–Gumukhola (920 m)	Dolakha	925	20740003	S. Noshiro et al.
9/13/2007	Dolakha (1750 m)–Nagdaha (890 m)–Malephu (940 m)–Gumukhola (920 m)	Dolakha	980	20740151	S. Noshiro et al.
9/13/2007	Dolakha (1750 m)–Nagdaha (890 m)–Malephu (940 m)–Gumukhola (920 m)	Dolakha	980	20740150	S. Noshiro et al.
9/13/2007	Dolakha (1750 m)–Nagdaha (890 m)–Malephu (940 m)–Gumukhola (920 m)	Dolakha	995	20740149	S. Noshiro et al.
9/13/2007	Dolakha (1750 m)–Nagdaha (890 m)–Malephu (940 m)–Gumukhola (920 m)	Dolakha	960	20740148	S. Noshiro et al.
9/13/2007	Dolakha (1750 m)–Nagdaha (890 m)–Malephu (940 m)–Gumukhola (920 m)	Dolakha	960	20740147	S. Noshiro et al.
9/13/2007	Dolakha (1750 m)–Nagdaha (890 m)–Malephu (940 m)–Gumukhola (920 m)	Dolakha	925	20740146	S. Noshiro et al.
9/13/2007	Dolakha (1750 m)–Nagdaha (890 m)–Malephu (940 m)–Gumukhola (920 m)	Dolakha		20740145	S. Noshiro et al.
9/13/2007	Dolakha (1750 m)–Nagdaha (890 m)–Malephu (940 m)–Gumukhola (920 m)	Dolakha	925	20740144	S. Noshiro et al.
9/13/2007	Dolakha (1750 m)–Nagdaha (890 m)–Malephu (940 m)–Gumukhola (920 m)	Dolakha	925	20740143	S. Noshiro et al.
9/15/2007	Jagat–Phedi–Buthu–Simigaun	Dolakha	1580	20740163	S. Noshiro et al.
9/15/2007	Jagat–Phedi–Buthu–Simigaun	Dolakha	1580	20740164	S. Noshiro et al.
9/15/2007	Jagat–Phedi–Buthu–Simigaun	Dolakha	1580	20740165	S. Noshiro et al.
9/15/2007	Jagat–Phedi–Buthu–Simigaun	Dolakha	1580	20740166	S. Noshiro et al.
9/15/2007	Jagat–Phedi–Buthu–Simigaun	Dolakha	1670	20740167	S. Noshiro et al.
9/15/2007	Jagat–Phedi–Buthu–Simigaun	Dolakha	1970	20740168	S. Noshiro et al.
9/22/2007	Dongang–Kelche–Suchek–Simigaun	Dolakha	2845	20740698	S. Noshiro et al.
9/16/2007	Simigaun (2030 m)–Suchek (kharka) (2490 m)–Kelche (2810 m)–Dongang (2780 m)	Dolakha	2825	20740381	S. Noshiro et al.
9/16/2007	Simigaun (2030 m)–Suchek (kharka) (2490 m)–Kelche (2810 m)–Dongang (2780 m)	Dolakha	2810	20740383	S. Noshiro et al.
9/16/2007	Simigaun (2030 m)–Suchek (kharka) (2490 m)–Kelche (2810 m)–Dongang (2780 m)	Dolakha	2810	20740384	S. Noshiro et al.
8/9/1999	Tamage-Banjan	Kaski	1950	9961182	M. Mikage et al.
8/1/1977	Shemma–Yakuwa–Lamobagar	Sankhuwasawa		770299	H. Ohashi et al.
8/2/1977	Lamobagar–Gola–above Shinbun	Sankhuwasawa		771940	H. Ohashi et al.
8/3/1977	Above Shinbun–HatiaGola	Sankhuwasawa		771970	H. Ohashi et al.
8/5/1977	HatiaGola–Digedanra–TaramBhanjyang–Honkon	Sankhuwasawa		773483	H. Ohashi et al.
9/3/1977	Sindua–Hile–Dhandkuta–TekuNala	Dhankuta		775560	H. Ohashi et al.
9/20/1963	Sundarijal Waterfall	Kathmandu	1600	6303366	H. Hara et al.
9/22/1963	Gokarna Forest	Kathmandu		6303367	H. Hara et al.
10/18/1963	Dhankuta	Dhankuta	1200	6303368	H. Hara et al.
10/19/1963	Dhankuta	Dhankuta	1200	6303369	H. Hara et al.
10/29/1963	Minchin Dhap–MulPokhari	Sankhuwasawa		6303370	H. Hara et al.
11/1/1963	Dumhan–Taplejung	Taplejung		6303371	H. Kanai et al.
11/10/1963	Helok–BaroyaKhimty	Taplejung		6303372	H. Hara et al.
7/7/1983	Suiket–Pathana	Kaski		8340051	H. Ohba et al.
7/5/1985	Shivalaya–Khasrubus–Deorali	Ramechhap		8540144	H. Ohba et al.
8/4/1985	Kyama–Patkare	Ramechhap		8571275	H. Ohba et al.
7/2/1988	Hile (Hille)–Pakhribas–GholiKharka (Dholikharka)	Dhankuta	1770	8840007	M. Suzuki et al.
7/2/1988	Hile (Hille)–Pakhribas–GholiKharka (Dholikharka)	Dhankuta	1830	8860007	M. Suzuki et al.
8/18/1988	Pokhara–Bagar–Hyangja–Suikhet–Majhbhatti–NaudandaPhedi (NaudharaPhedi)	Kaski	1160	8860430	M. Suzuki et al.
9/13/1988	Forest Office–Chhap–Shiwapuri Summit	Kathmandu	1620	8861045	M. Suzuki et al.
7/23/1990	Hile–Pakhribas–GholiKharka	Dhankuta	1800	9040023	M. Minaki et al.
8/13/1991	Baidep–Arun River–Num	Sankhuwasawa		9153461	H. Ohba et al.
5/11/1992	Mitlung–Thuma–Khokling	Taplejung	790	9261005	M. Suzuki et al.
5/14/1992	Papung–Dongen–Sewaden	Taplejung	2260	9261019	M. Suzuki et al.
5/24/1992	RamsyangPati–Jongim–WolangchungGola	Taplejung	2710	9261067	M. Suzuki et al.
6/5/1992	Amjilasa–bend–YanjolarchaKhola–Gyabla	Taplejung	2370	9261174	M. Suzuki et al.
6/6/1992	Gyabla–Nag PokariKhola–Fale–Ghunsa	Taplejung	2740	9261189	M. Suzuki et al.
6/14/1992	Yamphudin-SigreDanda	Taplejung	1900	9261278	M. Suzuki et al.
8/21/1998	Chyangrima–ChhokamKhola–Camp Site	Sankhuwasawa	2280	9840114	S. Noshiro et al.
9/23/1963	Balaju	Kathmandu			H. Hara et al.
7/26/2003	Godawari	Lalitpur	1570		T. Watanabe & K.J. Malla
11/24/2003	Sagakgaon	Gorkha			T. Watanabe & K.J. Malla

Barcode	Collection no.	Collection date	Collector(s)
E00240540	4441	9/16/1954	Stainton, J.D.A., W. Sykes and J. Williams
E00240656	1960	10/30/1954	Zimmermann, A.
E00293603	60/00	4/8/2000	W.R.Sykes
E00248246	A216	11/29/2004	Second Darwin Nepal Fieldwork Training Expedition
E00223626	9686196	9/21/1996	Botanical Expedition to West Nepal 1996
E00240549	1868	10/27/1954	Zimmermann A.
E00074569	9153461	8/13/1991	East Nepal Expedition to Koshi Zone (1991)
E00156584	9611796	1/19/1996	Botanical Expedition to Central Nepal, Chitwan (1996)
E00240539	3760	3/31/1952	Polunin, O., W. Sykes and J. Williams
E00240545	1868	10/27/1954	Zimmermann, A.
E00083102	9686196	9/21/1996	Mikage, M., R. Hirano, A. Takahashi and K. Yonekura
E00649206	280/99	4/27/1999	Sykes, W. R.
E00083140	9040023	7/23/1990	Botanical Expedition to Himalaya (1990)
E00250164	BY11	9/11/2005	Third Darwin Nepal Fieldwork Training Expedition

Locality	District	Altitude (m)
Aghor, Namtar-3	Makwanpur	1850
Bahundanda	Gorkha	1690
Chandragiri	Kathmandu	2050
Chitre	Parbat	2020
Dapcha	Kavre	1490
Dubichaur	Gulmi	1470
Godam	Kathmandu	1690
Gurba-3, Hastichaur	Gulmi	980
Harpukot	Gulmi	1120
Kalikot	Kalikot	2970
Godawari	Lalitpur	2230
Mahadevthan	Kathmandu	1580
Maldi-1, Hastichaur	Gulmi	1050
Naubise	Dhading	1270
Nayagaun, Daman-5	Makwanpur	2140
Panchase	Parbat	2830
Phidim	Phidim	1600
Phulchoki	Lalitpur	1760
Satghumti, Thakre-8	Dhading	1090
Thankot	Kathmandu	1950

Barcode	Collection no.	Collection date	Collector(s)	Cite as:(Stable URL)
E00240657	1960	10/30/1954	Zimmermann, A.	http://data.rbge.org.uk/herb/E00240657
E00240654	5243		Polunin, O., Sykes, W. & Williams, J.	http://data.rbge.org.uk/herb/E00240654
E00397283	B146	6/16/2008	Flora of Nepal Collecting Trip (2008)	http://data.rbge.org.uk/herb/E00397283

Pearson Correlations Coefficient (r)
Variables	bio_1	bio_2	bio_3	bio_4	bio_5	bio_6	bio_7	bio_8	bio_9	bio_10	bio_11	bio_12	bio_13	bio_14	bio_15	bio_16	bio_17	bio_18	bio_19	ai_ann	pet_ann
bio_1	1
bio_2	0.333**	1
bio_3	0.003	0.324**	1
bio_4	0.231	0.551**	−0.586**	1
bio_5	0.959**	0.459**	−0.143	0.460**	1
bio_6	0.980**	0.192	−0.014	0.121	0.927**	1
bio_7	0.342**	0.773**	−0.342**	0.936**	0.564**	0.214	1
bio_8	0.997**	0.364**	0.003	0.263*	0.960**	0.966**	0.373**	1
bio_9	0.988**	0.278*	0.013	0.175	0.938**	0.982**	0.280*	0.981**	1
bio_10	0.997**	0.370**	−0.039	0.302*	0.975**	0.970**	0.407**	0.997**	0.981**	1
bio_11	0.996**	0.288*	0.059	0.141	0.934**	0.987**	0.260*	0.990**	0.989**	0.986**	1
bio_12	0.402**	−0.119	0.235	−0.294*	0.253*	0.417**	−0.261*	0.396**	0.418**	0.362**	0.430**	1
bio_13	0.381**	−0.061	0.204	−0.218	0.247*	0.374**	−0.179	0.381**	0.389**	0.346**	0.401**	0.982**	1
bio_14	0.016	−0.252*	−0.253*	0.035	0.041	0.087	−0.085	0.011	0.017	0.026	0.018	−0.003	−0.107	1
bio_15	0.258*	0.217	0.114	0.086	0.189	0.146	0.17	0.283*	0.219	0.246*	0.243*	0.498**	0.628**	−0.502**	1
bio_16	0.396**	−0.061	0.225	−0.238	0.257*	0.392**	−0.194	0.393**	0.405**	0.358**	0.417**	0.989**	0.997**	−0.108	0.605**	1
bio_17	−0.07	−0.195	−0.621**	0.374**	0.093	0.009	0.222	−0.086	−0.061	−0.033	−0.099	−0.275*	−0.343**	0.733**	−0.614**	−0.344**	1
bio_18	0.222	−0.159	0.410**	−0.469**	0.041	0.238	−0.419**	0.216	0.244*	0.173	0.263*	0.944**	0.936**	−0.135	0.510**	0.943**	−0.448**	1
bio_19	−0.162	−0.106	−0.644**	0.476**	0.04	−0.099	0.321**	−0.176	−0.147	−0.117	−0.202	−0.374**	−0.420**	0.579**	−0.613**	−0.424**	0.961**	−0.528**	1
ai_ann	0.032	−0.407**	0.221	−0.505**	−0.144	0.074	−0.538**	0.025	0.057	−0.016	0.074	0.901**	0.882**	0.054	0.414**	0.881**	−0.252*	0.918**	−0.358**	1
pet_ann	0.967**	0.538**	0.055	0.362**	0.975**	0.922**	0.512**	0.970**	0.946**	0.976**	0.951**	0.302*	0.298*	−0.017	0.248*	0.309*	−0.05	0.123	−0.108	−0.096	1

*Correlation is significant at the 0.05 level (2-tailed).

**Correlation is significant at the 0.01 level (2-tailed).

GLC 2000 LULC Classification of Nepal
S.N.	GLC 2000 code	LULC type	Description
1	1	Tree cover, broadleaved evergreen, closed to open (>15%)	The main layer consists of broadleaved evergreen closed to open trees. The crown cover is between 100 and 15%. The height is in the range of >30–3 m but may be further defined into smaller range.
2	2	Tree Cover, broadleaved deciduous, closed (>40%)	The main layer consists of broadleaved deciduous closed to open trees. The crown cover is between 100 and 15%. The height is in the range of >30–3 m but may be further defined into smaller range.
3	4	Tree cover, needleleaved evergreen, closed to open (>15%)	The main layer consists of needleleaved evergreen closed to open trees. The crown cover is between 100 and 15%. The height is in the range of >30–3 m but may be further defined into smaller range.
4	11	Shrubcover, closed to open (>15%), evergreen (broadleaved or needle-leaved)	The main layer consists of broadleaved/needleleaved evergreen closed to open thicket. The crown cover is between 100 and 15%. The height is in the range of 5–0.3 m but may be further defined into smaller range.
5	12	Shrubcover, closed to open (>15%), deciduous (broadleaved)	The main layer consists of broadleaved deciduous closed to open thicket. The crown cover is between 100 and 15%. The height is in the range of 5–0.3 m but may be further defined into smaller range.
6	13	Herbaceous cover, closed to open (>15%)	The main layer consists of closed to open herbaceous vegetation. The crown cover is between 100 and 15%. The height is in the range of 3–0.03 m but may be further defined into smaller range.
7	14	Sparse Herbaceous or sparse Shrub cover	The main layer consists of sparse herbaceous vegetation. The crown cover is between (20–10) and 1%. The sparseness of the vegetation may be further specified. The main layer consists of sparse shrubs.
8	16	Cropland (upland crops or inundated/flooded crops as e.g. rice)	Primarily vegetated areas containing more than four percent vegetation during at least two months a year. The environment is influenced by the edaphic substratum. The vegetative cover is characterized by the removal of the (semi) natural vegetation and replacement with a vegetative cover resulting from human activities. This cover is artificial and requires maintenance. It is grown with the intention to be managed and/or (partly) harvested at the end of the growing season. Before or after harvest there may be a period without vegetative cover. The environment is significantly influenced by the presence of water over extensive period of time, i.e. water is present for more than three months a year and when water is present less than three months a year, it is present 75 percent of the flooding time.
9	19	Bare Areas	Primarily non-vegetated areas containing less than four percent vegetation during at least 10 months a year. The environment is influenced by the edaphic substratum. The cover is natural. Included are areas like bare rock and sands.
10	20	Water Bodies (natural or artificial)	The land cover consists of artificial water bodies. A further specification can be made in flowing or standing water. The land cover consists of natural water bodies.
11	21	Snow or Ice (natural or artificial)	The land cover consists of artificial snow. The land cover consists of artificial ice. A further specification can be made in moving or stationary ice. The landcover consists of snow and ice.
12	22	Urban Areas	The land cover consists of built up area(s).

Shapiro–Wilk Statistics (P-value)
Species		Model1	Model2	Model3	Model4
Alnus nitida	Training AUC	0.565	0.073	0.936	0.094
	Test AUC	0.008	0.053	0.614	1.28e−06
Alnus nepalensis	Training AUC	0.262	0.335	0.333	0.722
	Test AUC	0.488	0.483	0.389	0.248