Altitudinal Distribution and Habitat Suitability of 17 Inga Species (Fabaceae) in the Northern Ecuadorian Andes: A MaxEnt Modeling Approach for Conservation and Agroforestry Planning

1. Introduction

The genus Inga Mill. belongs to the family Fabaceae, subfamily Mimosoideae, and comprises approximately 300 species distributed from Mexico to southern Brazil and Bolivia [1,2]. This genus constitutes one of the most diverse and ecologically important in Neotropical forests, where its species perform essential functions such as biological nitrogen fixation, shade provision in agroforestry systems, and food production for wildlife [3,4]. Its pioneer character and tolerance of degraded soils also make it a widely used component of shade-grown agroforestry systems and land-restoration schemes across the Neotropics [3,5].

In Ecuador, the genus Inga exhibits notable diversity, with more than 60 reported species distributed from Amazonian lowlands to inter-Andean valleys [6,7]. Imbabura Province, located in the Northern Sierra of Ecuador (0°07–°52′ N, 77°45–°79°10′ W), presents exceptional environmental heterogeneity with altitudinal gradients ranging from 400 to 4939 m a.s.l., generating a diversity of microhabitats that favor the presence of multiple species of this genus [8,9].

Species Distribution Models (SDMs) are fundamental tools for understanding geographic distribution patterns and the environmental factors that determine species presence [10,11]. Among available algorithms, Maximum Entropy (MaxEnt) has demonstrated superior performance for modeling species distributions with presence-only data, particularly when the number of records is limited [12,13]. MaxEnt estimates the maximum entropy probability distribution subject to constraints derived from environmental conditions at known presence sites [14]. However, recent consensus guidelines for SDM practice emphasize the need for multi-metric model evaluation, variable selection to control multicollinearity, and explicit treatment of sampling bias to ensure model transferability and conservation relevance [15,16].

Despite the ecological and socioeconomic importance of Inga in the Ecuadorian Andes, studies on its detailed geographic distribution and the environmental factors determining its presence are scarce [17]. This lack of information limits the capacity to implement effective conservation and sustainable management strategies for these species [18] and hinders the development of evidence-based land-use planning that incorporates agroforestry as a tool for territorial management.

We hypothesize that (i) the altitudinal gradient and associated edaphoclimatic conditions are the primary determinants of Inga species distribution in Imbabura, and (ii) specialist species with restricted altitudinal ranges will yield higher model performance metrics (AUC, TSS) than generalist species, owing to their tighter environmental tolerances. The present study aims to: (i) document the geographic distribution of Inga species in Imbabura Province; (ii) generate and rigorously validate potential distribution models using MaxEnt; and (iii) identify the environmental variables of greatest influence on the distribution of each species, with the aim of contributing to conservation planning and agroforestry management in the region.

2. Materials and Methods

2.1. Study Area

Imbabura Province is located in the Northern Sierra of Ecuador, between coordinates 0°07′ and 0°52′ N latitude and 77°45′ and 79°10′ W longitude, covering an area of 4611 km². The study area encompasses six cantons (Ibarra, Otavalo, Cotacachi, Antonio Ante, Pimampiro, and Urcuquí) and presents an altitudinal gradient from 400 m in the subtropical zone of the biogeographic Chocó to 4939 m a.s.l. at Cotacachi volcano. This altitudinal variation generates five climatic belts: tropical, subtropical, temperate, cold, and páramo, with annual precipitation ranging from 250 to 3500 mm and mean temperatures between 6 °C and 24 °C [8,19].

2.2. Occurrence Data Collection and Spatial Thinning

Data on the occurrence of Inga species were obtained from two main sources: (a) a review of national and international herbarium databases, including the Herbarium of the Technical University of the North (HUTN), the National Herbarium of Ecuador at the National Institute of Biodiversity (QCNE), the Missouri Botanical Garden (MO), the Aarhus University Herbarium (AAU), and the Field Museum Herbarium (F) [20]; and (b) field expeditions conducted between 2024 and 2026, during which specimens were collected through targeted sampling in the province’s six cantons, under scientific collection permit No. MAE-DNB-CM-2024-0078 issued by the Ministry of Environment, Water, and Ecological Transition (MAATE). The geographic coordinates of all collected specimens were recorded using Kobotoobox. Duplicate records and those with imprecise coordinates (coordinate uncertainty > 1,000 m) were removed. In total, 181 georeferenced records of 17 Inga species were documented (Table S1)

To mitigate the spatial sampling bias inherent in herbarium records—which tend to cluster near roads and populated areas— To align with the 30 arc-second (~1 km²) resolution of the environmental predictors and avoid matching multiple occurrences to a single grid cell, we applied spatial thinning using the spThin package (v.0.2.0) in R [21], maintaining a minimum geographic distance of 1 km between retained records. This procedure reduced inter-point spatial autocorrelation while preserving maximum species coverage. After thinning, 10 species retained sufficient records (≥6 presence points) for MaxEnt modeling [22,23]. The complete georeferenced dataset is available through GBIF (https://doi.org/10.15468/dl.8sua8u).

2.3. Environmental Variables and Multicollinearity Screening

Five categories of environmental variables were initially compiled as predictors: (1) bioclimatic variables (Bio 1–19) from WorldClim 2.1 at 30 arc-sec resolution (≈1 km²) [24]; (2) topographic variables (elevation, Terrain Ruggedness Index TRI, slope) derived from the ALOS PALSAR digital elevation model at 12.5 m resolution [25]; (3) global soil variables (nitrogen, pH, texture, sand, silt, and clay content at 0–30 cm) from SoilGrids 2.0 at 250 m [26]; (4) national geopedological variables (soil order, geology) from the SIGTIERRAS program at 1:25,000 scale [27]; and (5) vegetation cover and productive systems from the Ministry of Environment, Water and Ecological Transition (MAATE) at 1:100,000 scale [28]. All layers were reprojected to WGS 84 coordinate system, resampled to a common 1 km² resolution using bilinear interpolation (continuous variables) or nearest-neighbor assignment (categorical variables) in ARCGIS pro, and clipped to the geographic extent of Imbabura Province

To address multicollinearity prior to modeling, Pearson correlation coefficients were calculated among all continuous predictor variables. Variables exhibiting |r| > 0.75 with another variable of higher ecological interpretability were removed from the predictor set [15]. Topographic variables (elevation, TRI, slope) and global soil attributes (SoilGrids) exhibited high multicollinearity (Pearson |r| > 0.75) with the dominant bioclimatic and geopedological predictors, and were consequently excluded to prevent overfitting. This screening retained a final set of 12 uncorrelated predictors per species (Table 1; Table S2), substantially reducing the risk of model overfitting. From the initial 19 bioclimatic variables, only Bio1, Bio3, Bio7, Bio8, Bio9, Bio12, and Bio17 were retained after screening; Tmax, Tmin, and water vapor pressure (VAPR) were also included as continuous predictors.

2.4. MaxEnt Modeling and Parameter Tuning

Potential distribution models were generated using the Maximum Entropy algorithm (MaxEnt v.3.4.4) [14,29]. For each species, regularization multiplier (RM) and feature class (FC) combinations were evaluated using the ENMeval package (v.2.0.4) in R [30], testing RM values of 0.5, 1.0, 1.5, 2.0, and 3.0 with FC combinations of L (linear), LQ (linear + quadratic), and LQH (linear + quadratic + hinge). The optimal configuration for each species was selected based on minimum corrected Akaike Information Criterion (AICc), favoring parsimony. For species with N ≤ 10, only L and LQ feature classes were considered to prevent overparameterization. Background was generated from 10,000 randomly distributed points within the province extent, with a maximum of 500 iterations and a convergence threshold of 0.00001. Logistic output format was used for all models. Final regularization multipliers and selected feature classes per species are reported in Table S3 (Supplementary Materials).

Model validation was performed using two complementary approaches. For species with N ≥ 10, data were partitioned into 75% training and 25% test sets with 10 bootstrap replicates; mean and standard deviation of AUC values are reported. For species with N < 10 (I. spectabilis, I. cocleensis, I. acuminata), leave-one-out (LOO) cross-validation was used to obtain a more robust estimate of model uncertainty [22]. Model performance was evaluated using: (1) the Area Under the ROC Curve (AUC); (2) the True Skill Statistic (TSS = Sensitivity + Specificity – 1) [31]; and (3) the Jackknife test of variable importance [32,33]. The TSS ranges from −1 (no better than random) to +1 (perfect model), and values > 0.6 are considered good for conservation applications [31]. Habitat suitability maps were classified into four categories: low (0–0.25), medium (0.25–0.50), high (0.50–0.75), and very high (0.75–1.0) [34].

Although three species in this study (I. spectabilis, I. cocleensis, and I. acuminata) presented a limited number of occurrence records (N = 7, 6–8, and 7, respectively), they were retained in the modeling process due to their high ecological and agroforestry relevance in the region. Model robustness for these species was ensured through LOO cross-validation, restriction to L and LQ feature classes, elevated regularization multipliers (RM ≥1.5), and evaluation using TSS in addition to AUC.

3. Results

3.1. Species Richness and Altitudinal Distribution

A total of 17 Inga species were documented in Imbabura Province, distributed along an altitudinal gradient spanning from 600 to 2800 m a.s.l. (Table 2). Of these, seven species were classified as habitat specialists, restricted to narrow altitudinal ranges (≤450 m amplitude), while three species were classified as generalists with altitudinal ranges exceeding 600 m; the remaining seven species exhibited intermediate distributions. The species with the highest number of records were I. insignis (28 records), I. densiflora (25), and I. punctata (20). Inga insignis presented the widest altitudinal distribution (1100–2800 m a.s.l.), while I. cinnamomea and I. silanchensis were restricted to low subtropical zones (~600 m a.s.l.). A complete list of georeferenced occurrence records per species is provided in Table S1 (Supplementary Materials).

3.2. MaxEnt Model Performance

MaxEnt models generated for the 10 species with sufficient records showed satisfactory to excellent predictive performance (Table 3). Training AUC values ranged from 0.9519 (I. cocleensis) to 0.9968 (I. punctata). Test AUC values ranged from 0.7727 (I. striata) to 0.9983 (I. feuillei). TSS values, calculated at the threshold maximizing the sum of sensitivity and specificity, ranged from 0.72 (I. striata) to 0.98 (I. feuillei), confirming adequate to excellent discrimination capacity across species. Six of the 10 modeled species had test AUC values above 0.95 and TSS values above 0.90. I. striata and I. insignis presented the lowest test AUC (0.7727 and 0.8795) and TSS (0.72 and 0.81) values, which can be attributed to their wide distribution and broad environmental tolerance as generalist species well-documented phenomenon in niche modeling literature [35,36]. For species modeled via LOO cross-validation (I. spectabilis, I. cocleensis, I. acuminata), mean AUC ± SD are reported across leave-one-out partitions.

3.3. Environmental Variable Contribution and Jackknife Evaluation

The variable contribution analysis and systematic Jackknife tests revealed highly differentiated, species-specific patterns among the modeled taxa (Table 4, Figure 1). The percentage contribution (%) and permutation importance (%) represent two complementary measures: the former reflects the incremental improvement during model training, while the latter evaluates the sensitivity of the model to random permutation of each variable in the test data and is therefore more robust for identifying ecologically meaningful predictors [32,33].

The systematic Jackknife evaluation (Table 4, Figure 1) highlighted distinct ecological limits. The distribution of I. edulis is strongly constrained by isothermality (bio_3), which achieved the highest predictive gain both when used alone and when omitted. Water vapor pressure (VAPR) was the most important single predictor for I. densiflora, while maximum temperature (Tmax) dominated models for I. feuillei. Temperature annual range (bio_7) constituted the dominant predictor for I. insignis and I. striata, consistent with their broad altitudinal distributions. Productive land-use systems and bio_17 (Precipitation of Driest Quarter) were critical for I. cocleensis, reinforcing a strong ecological association between this species and traditional agroforestry landscapes and mesic microhabitats. The strong response of agricultural land-use layers for I. edulis and I. cocleensis further reinforces their documented association with Andean agroforestry systems [3,5].

3.4. Potential Suitability Areas

Potential distribution maps revealed high suitability areas concentrated mainly in three zones of the province: (1) the Chota–Mira valley (600–1200 m a.s.l.) for subtropical species such as I. edulis, I. spectabilis, and I. cocleensis; (2) the mid-temperate belt (1400–2200 m a.s.l.) for I. densiflora, I. punctata, and I. venusta; and (3) high inter-Andean valleys (2200–2800 m a.s.l.) for I. insignis and I. striata. Areas of overlapping high suitability (≥40% of grid cells classified as high or very high) for multiple species simultaneously identified in the cantons of Cotacachi, Otavalo, and Ibarra, coinciding with traditional agroforestry system zones [37,38]. These overlap zones represent critical territories for integrated land-use planning that simultaneously promotes biodiversity conservation and sustainable agricultural production.

Figure 2. Potential habitat suitability maps for the 10 modeled Inga species in Imbabura Province, generated with MaxEnt v.3.4.4. Suitability categories: low (0–0.25), medium (0.25–0.50), high (0.50–0.75), and very high (0.75–1.0). White polygons indicate cantonal boundaries.

Figure 2. Potential habitat suitability maps for the 10 modeled Inga species in Imbabura Province, generated with MaxEnt v.3.4.4. Suitability categories: low (0–0.25), medium (0.25–0.50), high (0.50–0.75), and very high (0.75–1.0). White polygons indicate cantonal boundaries.

4. Discussion

The results of this study demonstrate that the distribution of Inga in Imbabura is strongly determined by the altitudinal gradient and associated edaphoclimatic conditions, consistent with patterns reported for other Fabaceae genera in the tropical Andes [39,40]. The identification of 17 species in a single province (4611 km²) represents notable richness, comparable to that reported for larger regions such as the Ecuadorian Amazon (>40 species in ~120,000 km²; [6,7]) and higher than that documented in equivalent altitudinal studies in the Colombian Andes (~12 species per province; [17]).

The high performance of MaxEnt models (AUC > 0.95 and TSS > 0.90 for 7 of 10 species) validates the utility of this algorithm for modeling tropical species distributions with limited presence data, as previously documented [12,13,22]. The combination of AUC and TSS as complementary metrics provides a more robust assessment than AUC alone: while AUC is a threshold-independent measure of overall discriminative capacity [16], TSS accounts simultaneously for commission and omission errors and is less sensitive to model prevalence, making it more suitable for conservation decision-making [31]. The relatively lower AUC and TSS values for I. striata (Test AUC = 0.7727; TSS = 0.72) and I. insignis (Test AUC = 0.8795; TSS = 0.81) reflect the inherent difficulty of modeling generalist species with wide environmental tolerance, a well-documented phenomenon in niche modeling literature [35,36] consistent with our initial hypothesis.

Temperature annual range (bio_7) emerged as the dominant predictor for I. insignis and I. striata, suggesting that these wide-ranging generalists respond primarily to broad thermal variation across their elevational gradients. In contrast, isothermality (bio_3) was most critical for I. edulis and I. spectabilis, indicating that thermal stability—rather than absolute temperature—determines their distribution. Maximum temperature (Tmax) dominated the I. feuillei model, consistent with its restriction to cooler inter-Andean belts above 2000 m a.s.l. Water vapor pressure (VAPR) was the primary predictor for I. densiflora and I. venusta, reflecting the importance of atmospheric humidity across their mid-elevation ranges.

The importance of productive land-use systems as a key predictor variable for several species (I. cocleensis: 39.35% contribution; I. densiflora: 20.70%) suggests a close association between Inga species and modified agricultural landscapes, reinforcing the historical role of this genus in traditional Andean agroforestry systems [3,5]. This relationship has direct implications for land-use planning, as it indicates that the persistence of several Inga species depends on the maintenance of traditional agroforestry practices that are being progressively abandoned [41,42].

Some methodological considerations should be acknowledged. First, three species were modeled with a restricted number of occurrence records ($\mathrm{N}=6-7$ for I. spectabilis, I. cocleensis, and I. acuminata). Although this sample size sits at the lower limit for spatial analysis, the integration of LOO cross-validation, stringent feature class constraints (L and LQ only), and conservative regularization multipliers ($\mathrm{R}\mathrm{M}\ge 1.5$) effectively controlled overfitting risk [22,30], rendering these models valuable preliminary tools for regional agroforestry zoning. Second, the spatial homogenization of predictors to a 1 km² grid—while necessary to match global bioclimatic baselines—unavoidably smooths fine-scale topographic and edaphic variations. Consequently, microhabitat dynamics for localized taxa like I. cinnamomea and I. silanchensis should be interpreted with caution until higher-resolution downscaled climatic surfaces become available for the region

The generated suitability maps identify priority areas for in situ conservation of Inga in Imbabura, particularly in altitudinal transition zones. These areas partially overlap with the Mira River Watershed Protective Forest and the Cotacachi–Cayapas Ecological Reserve [48,49]. Protection gaps were identified in the subtropical zones of the northwestern province, where I. cinnamomea and I. silanchensis present restricted distributions and may be vulnerable to ongoing deforestation [50,51]. Future studies incorporating Representative Concentration Pathways (RCP 4.5 and 8.5) for 2050 and 2070 would allow a quantitative assessment of range shifts and potential contraction zones for these Andean species.

5. Conclusions

This study documents the presence of 17 Inga species in Imbabura Province, Ecuador, providing the first detailed characterization of the geographic distribution of this genus in the northern Ecuadorian Andes. MaxEnt models, validated through multi-metric assessment (AUC, TSS, Jackknife) and complemented by multicollinearity screening (Pearson |r| < 0.75) and spatial thinning, demonstrated high predictive performance (Training AUC range: 0.9519–0.9968; TSS range: 0.72–0.98). Productive land-use systems, water vapor pressure, geopedological units, and bioclimatic temperature and precipitation variables (bio_3, bio_7, bio_12, bio_17) were the most influential factors. The classification of species into specialists and generalists—supported by both AUC and TSS evidence—provides a framework for prioritizing conservation actions focused on species with restricted distributions and high vulnerability. uture studies should incorporate climate change scenarios based on the Shared Socioeconomic Pathways (e.g., SSP2-4.5 and SSP5-8.5) from the CMIP6 framework, apply target-group background corrections, and explore connectivity modeling to identify priority corridors between fragmented suitable areas.