Identifying current and future species' distributions is a key element for planning management strategies and assessing conservation status in a rapidly changing environment (Miller 2010, Lawler et al. 2011). However, assessing geographic ranges of rare and understudied species can be problematic due to scarcity of available data, resulting in limited information for use by conservation managers (Pearce and Boyce 2006). Species Distribution Models (SDMs) can overcome deficiencies in distributional data by correlating the underlying environmental data from known occurrences and predicting to areas of highest environmental suitability based on the similarity to environmental conditions of known occurrence points (Scott et al. 2002, Pearce and Boyce 2006, Elith and Leathwick 2009). For example, SDMs can be used to direct survey efforts to find new populations of known species (Raxworthy et al. 2003), define protected areas and potential reintroduction sites (D'Elia et al. 2015, Bennett et al. 2017), and predict future effects of climate change (Aryal et al. 2016). Most traditional atlas range maps are too vague to direct action to specific areas within a species' distribution (Rodríguez et al. 2007). SDMs can provide a more focused approach for current and future climate scenarios, enabling targeted long-term management to core areas of highest environmental suitability where the focal species is most likely to occur currently and into the future.

Climate change is predicted to reduce the amount of suitable climate space for many narrow climate-adapted taxa, resulting in poleward latitudinal and upslope altitudinal range shifts (Walther et al. 2002, Parmesan and Yohe 2003, Chen et al. 2011). The predicted effects of increasing temperatures are most likely to affect bird species at higher latitudes or elevations, where poleward and upslope elevational range shifts with increased temperatures are predicted to be higher than at lower latitudes and elevations (Şekercioğlu et al. 2008, Sodhi et al. 2011, Freeman et al. 2018). Conversely, within the tropics, where precipitation may be the most important determinant of avian distributions (Pearce-Higgins and Green 2014), climate change may not necessarily result in poleward range shifts but in multidirectional shifts in suitable climate space (VanDerWal et al. 2013).

The predicted effects of climate change are thus unclear for species such as the Laggar Falcon (Falco jugger), which occupies a restricted range of arid climatic conditions within tropical and subtropical open habitats such as deserts, steppes, and savannas (Naoroji 2006). This species is a medium-sized raptor within the sub-genus Hierofalco, along with other congeneric “desert falcons”: Gyrfalcon (Falco rusticolus), Lanner Falcon (Falco biarmicus), and Saker Falcon (Falco cherrug; Nittinger et al. 2007). The Laggar Falcon has an Indo-Malayan distribution from southeastern Iran and Afghanistan in the west, across the Indian subcontinent to Myanmar in the east (Ferguson-Lees and Christie 2005, Naoroji 2006). Therefore, the distribution of the Laggar Falcon straddles the Tropic of Capricorn and extends into the Tropical zone across lower elevations. This falcon is regarded as a climatic semi-specialist, adapted to arid and semi-arid climates (Naoroji 2006, Finlayson 2011), from sea level up to 1000 masl in open, sparsely wooded landscapes with an estimated Extent of Occurrence (EOO) covering 4.2 million km² (Ferguson-Lees and Christie 2005, BirdLife International 2016). However, current distribution maps lack detail (BirdLife International 2016), resulting in a significant knowledge gap for defining current distribution and extrapolating to future areas of highest climatic suitability.

The Laggar Falcon is the least known of all the Hierofalcons (Cade 1982). Indeed, a recent review found that the species has been the subject of only four peer-reviewed studies, compared to 102 for Gyrfalcon, 70 for Saker Falcon, and 41 for Lanner Falcon (Buechley et al. 2019). Currently, Laggar Falcon populations are thought to be declining rapidly, and the species is categorized as Near Threatened on the IUCN Red List, due to increased pesticide use and incidental trapping for the falconry trade (BirdLife International 2016). Although not straightforward, predicting where future suitable climate space will exist for Laggar Falcons might assist long-term conservation management for this species.

With reliable distribution estimates based on climatic variables, the effect of climate change on Laggar Falcon populations can be assessed, and conservation actions implemented to assist in species management. Further, an initial assessment of climatic factors influencing Laggar Falcon distribution will highlight climatic constraints, and provide fundamental biogeographical information for the species. Our study aimed to increase understanding of the biogeographical constraints on distribution of the Laggar Falcon for use within conservation actions and management. The specific research objectives were to: (1) update current distribution for the Laggar Falcon and define core areas of climatic suitability, (2) determine the key climatic constraints that explain these distributions, and (3) predict future distribution under various climate change scenarios. Specifically, we expect range size of the Laggar Falcon to contract and shift multidirectionally in line with future climate change projections.

Methods

Laggar Falcon Occurrence Data. Presence-only occurrence data for the Laggar Falcon were sourced from the Global Raptor Impact Network (GRIN; The Peregrine Fund 2018), an information system containing data for all raptor species. For Laggar Falcons, GRIN contains occurrence data from the Global Biodiversity Information Facility (2019), which mostly consist of eBird records (89.5%, Sullivan et al. 2009). GRIN also contains three records of Laggar Falcon nests collected by the authors (SK and RS). We cleaned occurrence data by removing duplicate records and those with no georeferenced location. To account for spatial autocorrelation, sampling bias, and clustering in the occurrence points, we selected a spatial filter of 10 km from each occurrence point to minimize the effects of over-sampling in highly surveyed areas, using the “thin” function in the R package SPTHIN (Aiello-Lammens et al. 2015). The relatively coarse thinning distance of 10 km was chosen based on the environmental heterogeneity of the study area, and high potential of sampling bias from citizen science datasets (Beck et al. 2013). Spatial filtering by removing clustered occurrence points reduces model over-fitting, minimizes omission errors, and improves model predictive performance (Kramer-Schadt et al. 2013, Boria et al. 2014, Radosavljevic and Anderson 2014), and performs better than other spatial bias correction methods (Fourcade et al. 2014).

We measured spatial autocorrelation in occurrence data using Global Moran's I index on an inverse Euclidean distance matrix projected into Albers Equal Area for India. Moran's I is an index ranging from –1 to +1, with values closer to zero indicating no spatial autocorrelation, negative values indicate negative spatial autocorrelation, and positive values indicating positive spatial autocorrelation. We measured spatial clustering using Nearest Neighbor Index (NNI) in the R package SPATIALECO (Evans 2017) with a convex hull window. NNI is the ratio of the observed distance divided by the expected distance between neighbors in a hypothetical random distribution. NNI values <1 indicate spatial clustering, values >1 indicate random dispersion, and those closer to 2 indicate regular dispersion.

Climatic Predictors. Bioclimatic data for current distributions were obtained from WorldClim (version 1.4), generated through interpolation of average monthly weather station climate data from 1960 to 1990 (Hijmans et al. 2005). We used solely climatic predictors, rather than other topographical or habitat variables, as an indication of changing climatic conditions assuming all else remains equal. Although there are many factors that affect species' distributions, reliable SDMs for future changes in distribution can only be developed based on future climate model predictions. At present we are unable to predict how habitat will change by 2050 but there are a suite of future global climate change models that can be used to make predictions based on current climate constraints. Further, at broad scales examined in this study, climate is viewed as the main driver of species' distributions and thus bioclimatic predictors are the most appropriate variables to use (Pearson and Dawson 2003). We cropped raster layers ≤7° over the extent of Laggar Falcon occurrences using a delimited polygon. This improves model predictive power by defining the full range of suitable environmental conditions across the known range of the Laggar Falcon (Lawler et al. 2011) and reducing the background area used for testing points used in model evaluation (Radosavljevic and Anderson 2014).

Multicollinearity between environmental predictor variables can bias models by overrepresenting the biological relevance of correlated variables (Franklin 2009). Before model construction, we tested all 19 bioclimatic variables for multicollinearity using Variance Inflation Factor (VIF) analysis (Guisan et al. 2006, Hair et al. 2006) in the R package USDM (Naimi et al. 2014). VIF is based on the square of multiple correlation coefficients, regressing a single predictor variable against all other predictors. VIF tests can detect hidden correlations in predictors not always apparent in pair-wise correlations. VIF >10 indicates collinearity in the variables; thus, we used a stepwise elimination of highly correlated variables, retaining predictors with a more stringent VIF threshold of <5 (see   Table S1 (1_rapt-54-01-01_s01.pdf) in Supplemental Materials, available online), considered as suitable for multivariable correlation (Dormann et al. 2013). We then checked the remaining variables for collinearity using Pearson's Correlation Coefficient, with only variables r = |0.7| retained for consideration as predictors. After removing highly correlated variables, we included eight bioclimatic variables as predictors (mean diurnal temperature range, isothermality, mean temperature wettest quarter, mean temperature driest quarter, precipitation driest month, precipitation seasonality, precipitation warmest quarter, precipitation coldest quarter), at a spatial resolution of 2.5 arc-min (approximately 4.5 km² resolution at the equator).

We based final predictor selection on representing seasonal climatic trends, extremes, and limiting environmental factors strongly related theoretically and empirically to species' distributions (Stockwell 2006, Reineking et al. 2016, Bradie and Leung 2017), and specifically to distributions of vagile bird species in arid environments (Reside et al. 2010). Because current knowledge of Laggar Falcon biology is limited, we selected predictors known to affect avian distributions in arid environments (Sutton and Puschendorf 2018). We selected mean diurnal temperature range as an important predictor for diurnal species in arid environments because variability in daily temperatures likely affects avian survival and population viability (Briga and Verhulst 2015). In arid and semi-arid ecosystems, rainfall regime and seasonality is predicted to have a strong effect on avian survival, food availability, and reproductive effort (Winterbottom and Rowan 1962, Dean et al. 2009). We thus selected predictors reflecting seasonal and monthly extremes of precipitation interacting with temperature as potential limiting factors on Laggar Falcon distribution based on existing knowledge from other arid environment avian species (Gargett et al. 1995, Zann et al. 1995, Illera and Diaz 2006, Dean et al. 2009, Cavalcanti et al. 2016).

For future predictions, we used four General Circulation Models (GCMs) from the Coupled Model Inter-comparison Project Phase 5 (CMIP5) to predict Laggar Falcon distributions for two future climate scenarios in 2050 (averaged over the time period 2041–2060): INMCM4 (Institute of Numerical Mathematics Climate Model, v4.0, Volodin et al. 2010), MIROC5 (Model for Interdisciplinary Research on Climate, v5, Watanabe et al. 2010), MPI-ESM-LR (Max Planck Institute Earth System Model Low Resolution, Giorgetta et al. 2013), and NorESM1-M (Norwegian Earth System Model, Bentsen et al. 2013). We selected four GCMs to account for variation in model outputs, and any uncertainty in single model predictions (Lutz et al. 2016), with all four GCMs predicting future climate for south Asia well (Mishra et al. 2014, Sharmila et al. 2015). We aimed to select a representative choice of the least interdependent climate models spanning the spectrum of variability in future climate change prediction (Sanderson et al. 2015). Data were downloaded from the WorldClim database (version 1.4, Hijmans et al. 2005) for two different CMIP5 greenhouse gas concentration scenarios or Representative Concentration Pathways (RCPs, Meinshausen et al. 2011): RCP4.5 (lower emissions), and RCP8.5 (high emissions). RCP 4.5 represents CO₂ concentrations stabilizing to 650 ppm by 2100 without exceeding that value (Thomson et al. 2011), and RCP 8.5 corresponds to CO₂ concentrations of 993 ppm to 2150 (Meinshausen et al. 2011) or “business as usual” (Riahi et al. 2011).

Species Distribution Models. We fitted SDMs using a maximum entropy machine-learning algorithm, MAXENT (version 3.3.3k, Phillips et al. 2006). MAXENT uses presence-background data, and is a high-performing SDM algorithm, more accurate than other SDM methods, even at low sample sizes (Gibson et al. 2007, Duan et al. 2013, Elith et al. 2006, Wisz et al. 2008). MAXENT correlates environmental variables underlying species occurrences against a random sample of background environmental conditions. Over-fitting in complex models is reduced by using L1-regularization (β, Hastie et al. 2005), further reducing bias in correlated predictor variables (Phillips and Dudík 2008, Elith et al. 2011). Within the MAXENT software, we selected logistic output as a continuous index of climate suitability, with 0 = low suitability and 1 = high climatic suitability. We used default model parameters for generated background absences (10,000) and convergent threshold (10^–5), and increased iterations to 5000 from the default 500 allowing for model convergence.

Optimal-model selection was based on Akaike's Information Criterion (Akaike 1974) corrected for small sample sizes (AIC_c, Hurvich and Tsai 1989) to determine the most parsimonious model from two key MAXENT parameters: regularization multiplier and feature classes (Warren and Seifert 2011). Tuning MAXENT parameters smooths response curves, limits sampling bias, and reduces over-fitting in presence-only predictions (Merow et al. 2013, Radosavljevic and Anderson 2014). Forty-eight candidate models of varying complexity were built by comparing a range of regularization multipliers from 0.5 to 4.0 in 0.5 increments, and five feature classes (Linear, Quadratic, Hinge, Product, Threshold) in all possible combinations using the “block” method of cross-validation (k = 5) within the ENMEVAL package in R (Muscarella et al. 2014). Block partitioning masks the geographical structure of the data according to latitude and longitude lines, dividing all occurrences into four spatially independent bins of equal numbers. By masking the geographical structure of test-data, the models are projected onto an evaluation region not included in the calibration process. All occurrence and background test points are assigned to their respective bins dependent on location, thus further reducing spatial autocorrelation between testing and training localities (Muscarella et al. 2014, Radosavljevic and Anderson 2014). We chose the block method as our preferred data-partitioning method, because we are transferring model predictions in time where there is a high possibility of encountering no-analog climate conditions.

Model Evaluation. We evaluated model performance within ENMEVAL using both threshold-independent and threshold-dependent measures (Radosavljevic and Anderson 2014). Area Under the Curve (AUC) of the Receiver Operating Characteristic plot (ROC) (Hanley and McNeil 1982), is a nonparametric, threshold-independent measure representing an overall value of model performance across all thresholds (Peterson et al. 2011), with AUC = 1.0 being the maximum predictive performance, and an AUC = 0.5 being no better than a random prediction (Fielding and Bell 1997). AUC values >0.9 are considered of high predictive accuracy (Franklin 2009), though a moderate AUC value >0.7 is deemed viable for conservation planning (Pearce and Ferrier 2000, de Carvalho et al. 2017). Further, to quantify any over-fitting present in the models, AUC_DIFF (AUC_TRAIN –AUC_TEST) was also used to measure predictive performance (Muscarella et al. 2014). This value is expected to be close to zero in low over-fitted models (Warren and Seifert 2011). Two threshold-dependent measures were employed based on omission rates from two threshold rules: minimum training presence (MTP) and 10% training presence (10%TP) omission rate thresholds. Omission rates report the proportion of training points that are outside of the model when converted into a binary prediction. Omission rates evaluate discriminatory ability and over-fitting at specified thresholds. Overall lower omission rates show improved discrimination between suitable and unsuitable areas (indicating higher performance), while over-fitted models show higher omission rates than expected by theory (Radosavljevic and Anderson 2014). Therefore, for models with low over-fitting the expectation in MTP is a value close to zero and for 10%TP a value close to 0.10.

AUC has been criticized as a measure of model performance for presence-background SDMs (Lobo et al. 2008, Jiménez-Valverde 2012). Therefore, we tested our final model prediction against random expectations using partial ROC (pROC), which estimates model performance by giving precedence to omission errors over commission errors (Peterson et al. 2008). We set function parameters with a 10% omission error rate, and 1000 bootstrap replicates on 50% test data to determine significant (α = 0.05) pROC values >1.0 in the R package ENMGADGETS (Barve and Barve 2013). Further, we used the null model approach from Raes and ter Steege (2007), which calculates 95% AUC Confidence Intervals (CI) on a frequency histogram, randomly drawing points without replacement on 999 replicate models under the same environmental parameters used in our best-fit model. With fitted 95% CI AUC values, we assessed whether our model accuracy was significantly higher than expected by chance (P < 0.05). We used continuous Boyce index (B) as a second threshold-independent evaluation metric on the final model prediction (Hirzel et al. 2006, Ramírez-Albores et al. 2016). B measures how much climate suitability predictions differ from a random distribution of observed presences across the predictive maps (Boyce et al. 2002). It is consistent with a Spearman correlation (r_s) with values of B ranging from –1 to +1, with positive values indicting climatic suitability predictions consistent with observed presences, and values closer to zero no different than a random model, and negative values indicating areas with frequent presences having low climatic suitability. B evaluation was used on all presence data points split into 80% training (B_TRAIN), 20% testing (B_TEST), plus a measure of over-fitting in the final models (B_DIFF). B was calculated using the default settings in the ECO-SPAT package in R (Di Cola et al. 2017), with a moving window for threshold-independence and 10 defined bins. Final continuous predictive maps used all the occurrence data points to achieve the highest predictive accuracy for estimating distribution (Fielding and Bell 1997, Fehérvári et al. 2012). Lastly, we used a jack-knife test to estimate variable performance within the optimal calibration model by excluding each value, then developing the model with a sole variable to determine percentage contribution and regularized training gain of each environmental variable to model performance.

For future predictive models, we calculated pairwise niche overlap metrics for all continuous MAXENT predictions to quantify how predictions from the four GCMs differ in geographic space using Schoener's D (Schoener 1968, Warren et al. 2008), and the I similarity statistic (Warren et al. 2008). Both statistics range from 0 (no overlap) to 1 (identical predictions), with Schoener's D focusing on niche similarity between habitats, whereas the I similarity statistic uses Hellinger distances that carry no biological assumptions but solely compares the probability distributions. Future continuous predictions were reclassified as binary presence/absence predictions to calculate future climate space and aid interpretation, with a 10% TP threshold to determine the number of incorrect occurrences with environmental conditions. We chose 10%TP as a suitable threshold for a small sample of presence-only occurrences to reject the lowest 10% of predicted values accounting for any uncertainty in the occurrence data (Pearson et al. 2007). Further, using 10%TP in MAXENT models generally maintains a higher proportion of presences correctly predicted by artificially reducing sample size (Pearson et al. 2007), with improved predictive ability if a 10%TP threshold is set. Model development, GIS analysis, and predictive maps were built in R (version 3.5.1, R Core Team 2018) using the DISMO (Hijmans et al. 2017) and RASTER (Hijmans 2016) packages.

Results

Laggar Falcon Occurrence Data. A total of 171 Laggar Falcon geo-referenced records were sourced from the GRIN database for inclusion in model calibration. Cleaned occurrence data were spatially autocorrelated (Moran's I = 0.750, P ≤ 0.001). After applying the 10-km spatial filter, spatial autocorrelation was reduced (Moran's I = 0.504, P ≤ 0.001), resulting in 115 filtered occurrence records for use in the calibration models. Cleaned occurrence data showed spatial clustering (NNI = 0.482, z = –12.969, P ≤ 0.001), with clustering moving towards random dispersion after applying the 10-km spatial filter (NNI = 0.630, z = –7.587, P ≤ 0.001).

Species Distribution Models. We selected the MAXENT feature class parameters Linear, Quadratic, Hinge, and a regularization multiplier of β = 2 based on the optimal-model output from all 48 candidate models (ΔAIC_c = 0.0, w_i = 0.673). All evaluation metrics from the block five-fold cross-validation mean test data had high predictive performance (AUC_TRAIN = 0.882, AUC_TEST = 0.802), with low over-fitting and good discrimination ability (AUC_DIFF = 0.080, MTP = 0.089, 10%TP = 0.220). Our best-fit model had significantly higher accuracy than expected by chance, placed within the highest 5% against 999 null models (P < 0.05). Testing our final predictive model against a random expectation resulted in a high mean pROC value of 1.716 (±0.082, range 1.432–1.929, P ≤ 0.01). Boyce index values in the final model prediction showed high positive correlation between predicted climatic suitability and presence training (B_TRAIN = 0.963) and test occurrence points (B_TEST = 0.911), with low over-fitting (B_DIFF = 0.052), demonstrating well-calibrated models.

Current Distribution. The current continuous predictive map using all occurrence data was consistent with the known distribution of the Laggar Falcon, but also predicted new areas with high climate suitability, and a core range centered in west-central India and southeastern Pakistan (Fig. 1,  S1 (1_rapt-54-01-01_s01.pdf) in Supplemental Materials, available online). In India, Gujarat and Rajasthan had highest climatic suitability, including most of the Thar desert and extending into western Madhya Pradesh. Further high climate suitability was also indicated along a narrow coastal strip of the Western Ghats from Gujarat, through the provinces of Maharashtra and Goa and into northern Karnataka. We also identified small pockets of high climate suitability in southern India outside of the current known range for the Laggar Falcon: the Kurnool and Mahbubnagar regions in west-central Andhra Pradesh (Fig. 1). In Pakistan, highest climate suitability was identified in the southern half of Sind province, bordering similar areas of high climate suitability in Gujarat and Rajasthan in India.

Figure 1.

Predicted distribution model for the Laggar Falcon. Final predictive map using all occurrence data for increased accuracy. Map denotes continuous logistic prediction with hotter red areas (values closer to 1) having higher climatic suitability. White filled circles indicate known Laggar Falcon occurrences. Continuous Boyce index (B) using 20% test data: B_TEST = 0.911.

Environmental Predictors. The distribution of the Laggar Falcon was most associated with mean temperature of the wettest quarter (BIO8), which contributed the highest percentage to model prediction (50.1%, Table 1), followed by precipitation seasonality (BIO15, 17.6%), and precipitation during the driest month (BIO14, 12. 2%, Table 1). The remaining climate predictors all contributed <6% to model prediction. Response curves to climate variables demonstrated pronounced thresholds for Laggar Falcon climate suitability (Fig. 2). Mean temperature in the wettest quarter plateaued at a climate suitability of 30–40°C, with suitable conditions of mean temperature in the driest quarter ranging from 10–15°C. Precipitation in both the warmest and coldest quarters, and in the driest month had highest climate suitability at zero rainfall, though Laggar Falcons can tolerate a higher level of seasonal precipitation of between 15–25 mm. Mean diurnal temperature range was high, peaking at 15°C, reflecting the mainly arid environments the falcons inhabit, with high variation in day-night temperatures. Jack-knife tests of variable importance showed that mean temperature in the wettest quarter, and precipitation seasonality had the highest regularized training gain (Fig. 3). Mean temperature in the wettest quarter had the highest gain when used in isolation, and so had the most useful information on suitable climatic conditions when used alone. Isothermality had the lowest decrease in gain when omitted, and therefore has the most information that is not present in the other variables to explain the climatic niche of the Laggar Falcon.

Table 1.

Percent contribution for bioclimatic variables used as environmental predictors in species distribution model for the Laggar Falcon. All temperature variables in °C and precipitation variables in mm.

Figure 2.

Response curves for each environmental variable model used as predictors in distribution model for the Laggar Falcon.

Figure 3.

Variation in regularized training gain for climatic predictors from a species distribution model for the Laggar Falcon using a jack-knife test of variable importance. Light gray bars represent regularized gain without the variable and dark gray bars represent regularized gain with only the one variable. All temperature variables (°C) and precipitation variables (mm).

Future Distributions. Predicted future distributions varied spatially among GCMs in where suitable climate space will exist in 2050 for Laggar Falcons, with wide variation among GCMs in estimated gain and loss of future suitable climate space. All future Laggar Falcon distributions from the four GCMs showed multidirectional shifts in both RCP scenarios, especially in the high emissions scenario, but with highest climate suitability still centered in Gujarat and Rajasthan in India and Sind, Pakistan (Fig. 4, 5). Niche overlap metrics showed high geographic overlap among all GCM predictions in both emission scenarios (Table 2). Future binary distributions for 2050 based on the 10% training presence threshold rule and two different climate scenarios averaged over the four GCMs, showed a mean gain of +6.0% (SE = 8.0%) to 2,096,160 km² (SE =158,513.4 km²) of suitable climate space in the RCP 4.5 scenario, and a mean reduction of 5.0% (SE = 8.5%) to 1,880,017 km² (SE = 167,620.4 km²) of suitable climate space in the RCP 8.5 scenario (Table 3,  Fig. S2, S3 (1_rapt-54-01-01_s01.pdf) in Supplemental Materials, available online), compared to the current predicted suitable climate space of 1,978,070 km².

Figure 4.

Predicted distribution in 2050 for the Laggar Falcon using four General Circulation Models (GCMs) under a low carbon emissions scenario (RCP4.5). Maps denote continuous logistic prediction with hotter red areas (values closer to 1) having higher climatic suitability.

Figure 5.

Predicted distribution in 2050 for the Laggar Falcon using four General Circulation Models (GCMs) under a high carbon emissions scenario (RCP8.5). Maps denote continuous logistic prediction with hotter red areas (values closer to 1) having higher climatic suitability.

Discussion

The geographic range of a species is a fundamental unit within biogeography (Brown et al. 1996), with SDMs having significant potential to make key contributions to spatial conservation planning (Lawler et al. 2011, Guisan et al. 2013). Our predictive distribution models for the Laggar Falcon demonstrated a core range in western India and southeastern Pakistan, consistent with an abundant center distribution (Andrewartha and Birch 1954, Brown et al. 1995), though use of only climatic variables may oversimplify complex biogeographical patterns (Dallas et al. 2017). Climatic suitability for the Laggar Falcon is highest in the center of its range, becoming less suitable towards range edges, defining priority areas for conservation management actions. Model predictive performance was high using a range of evaluation metrics and explained fundamental information on the climatic constraints on Laggar Falcon distribution across the Indian subcontinent. Predicted future distributions under two climate change scenarios demonstrated multidirectional range changes in suitable climatic space of up to 6% expansion in a lower carbon emissions scenario and up to 5% contraction for a high carbon emissions scenario. The SDMs given here for a rare and poorly understood species accurately predict current and estimate future areas of highest climatic suitability where management actions should be given precedence.

Revised Current Distribution. Our current distribution model updates previous distribution maps for the Laggar Falcon (Ferguson-Lees and Christie 2005, Naoroji 2006, Nittinger 2007, BirdLife International 2016), and captures the spatial complexity in the distribution of this rare raptor. According to our analysis, Gujarat and Rajasthan in India and southern Sind in Pakistan are the three regions containing the most suitable climatic conditions. We recommend Gujurat, Rajasthan, and Sind as high priority areas for population surveys and regular monitoring, but also recommend exploratory surveys into the Western Ghats and west-central Andhra Pradesh. At present there are no known conservation actions underway for the Laggar Falcon (BirdLife International 2016); thus, based on these models we recommend establishing surveys and monitoring Laggar Falcon populations in the priority areas described above. We focused solely on climatic conditions but recognise that our models have limitations. Distribution model predictions at this scale would likely be improved by including topographical predictors (Meineri and Hylander 2017), along with biotic interactions with potential competitors such as the Red-necked Falcon (Falco chicquera), distribution of prey resources (Wisz et al. 2013), and human pressures (Engler et al. 2017). Future modelling efforts should examine these other possible drivers of Laggar Falcon distribution. However, diagnostics revealed that our models perform well, despite omission of other potentially important variables.

Climatic Constraints. Defining climatic variables that drive distributions is important in the context of designing management plans that adapt to current and future climate scenarios (Prato 2012). Our model response curves and jack-knife tests clearly identify those key climatic characteristics that define the distribution of Laggar Falcons. This species can tolerate high mean temperatures above 30°C during the wettest quarter, but this tolerance decreases to a mean of 10°C during the driest quarter, suggesting a link between optimal seasonal temperature and precipitation. Though Laggar Falcons can withstand variable levels of seasonality in precipitation up to 25 mm, they generally prefer areas with zero or low rainfall in seasonal or monthly predictions, as expected from a species adapted to a semi-arid environment. Laggar Falcons have a broad tolerance of up to 15°C in both mean diurnal temperature range and isothermality, which points toward distribution being influenced by large temperature fluctuations on a daily and monthly basis relative to the year.

Table 2.

Niche overlap test metrics for continuous prediction maps for the Laggar Falcon using low (RCP4.5) and high (RCP8.5) emissions climate change scenarios from four General Circulation Models (GCMs).

Table 3.

Change in suitable future climate space for the Laggar Falcon using low (RCP4.5) and high (RCP8.5) emissions climate change scenarios from four General Circulation Models (GCMs).

Future Distributions. Climate change is expected to alter distributions of many bird species globally (Huntley et al. 2006), though responses of individual species will likely vary dependent on current status, location, and species' ecology (Møller et al. 2006). As expected, our future distribution models predict multidirectional range shifts in suitable climate space. Unexpectedly, suitable climatic conditions may expand 6% by 2050, but with a 5% contraction by 2070. However, there was high variation in the predicted amount of suitable future climate space among GCMs, though all were consistent in maintaining future suitable climate space in areas similar to the current core range. Laggar Falcons currently face multiple threats, although our models indicate climate change may not significantly affect populations. However, because of the small sample size used for our analysis, we urge caution in dismissing the threat of climate change for this species. Our results demonstrate how individual species may respond differently to changing climate, and that species-specific SDMs need to be calibrated following current best practice to effectively predict future distributions (Hijmans and Graham 2006, Elith et al. 2010). Using a range of future climate projections from multiple GCMs can account for any variation in future predictions, providing more realistic projections of future range shifts (Sanderson et al. 2015).

Recommendations for Monitoring. By identifying areas of optimal climate, monitoring and conservation efforts can be focused to track changes in population levels associated with environmental and climate change (Brown et al. 1995). However, we recognize that our analyses are based on relatively few occurrence samples, and thus may not capture the full suite of suitable environmental conditions (Sagarin and Gaines 2002). Therefore, where possible we recommend future sampling from across the entire geographic range of the Laggar Falcon, along with exploratory surveys into areas of highest environmental suitability. These data can then be fed back into our models, enabling more accurate range maps, and quantification of environmental constraints on distribution. Increasing the number of Laggar Falcon occurrence records within GRIN and establishing a network of field observers would be an initial first step in this process. With enhanced model predictions and niche quantification, we would have increased confidence in conservation actions such as recommending protected areas.

Climate envelope models have developed rapidly over the past 10 years and have shown their value as analytical tools in conservation biology (Pearson and Dawson 2003, Araújo and Peterson 2012). Our models accurately predict the current distribution of Laggar Falcons while also identifying climatic constraints and predicting future range shifts. We demonstrate that using online biodiversity data within a presence-only SDM framework can predict current distribution accurately and identify key biogeographical constraints to distribution. We further show that projecting current predictions into future climate change scenarios can yield useful results when applied in a robust modelling methodology. Globally, the most pressing research need for raptors is determining population trends (McClure et al. 2018), and the Laggar Falcon is no exception. This study should inform future monitoring efforts by refining the known distribution of the species and identifying areas predicted to experience changes in population levels within the next few decades.

Supplemental Materials.  Table S1 (1_rapt-54-01-01_s01.pdf). Selected bioclimatic variables used as predictors with VIF scores.

 Figure S1. (1_rapt-54-01-01_s01.pdf) Projected distribution into core area of highest suitability in Gujarat and Rajasthan, India, and Sind, Pakistan.

 Figure S2. (1_rapt-54-01-01_s01.pdf) .Reclassified binary predictive maps using 10% training threshold under a low carbon emissions scenario (RCP 4.5).

 Figure S3. (1_rapt-54-01-01_s01.pdf) Reclassified binary predictive maps using 10% training threshold under a high carbon emissions scenario (RCP 8.5).