Combining occurrence and abundance distribution models for the conservation of the Great Bustard

Introduction

The knowledge of species occurrence and abundance distribution provides fundamental information for conservation biology (VanDerWal et al., 2009; Drew, Wiersma & Huettmann, 2011; Primack, 2012; Johnston et al., 2015). Understanding how environmental factors are related to species occurrence and abundance distribution and how they are explicit in time and space are priorities in current biodiversity conservation (Drew, Wiersma & Huettmann, 2011; Martín et al., 2012).

Species distribution models (SDMs) are empirical ecological models that relate species observations to environmental predictors (Guisan & Zimmermann, 2000); usually this process is done using machine learning algorithms (Drew, Wiersma & Huettmann, 2011, see Mi et al., 2017 for an application). SDMs have become important and essential tools in ecology, biogeography, climate change research, conservation, and management because of on their spatial occurrence prediction capacities (Peterson et al., 2002; Guisan & Thuiller, 2005; Elith et al., 2006; Araújo & New, 2007; Mi, Huettmann & Guo, 2016). SDMs built with count data are called species abundance models (SAMs) (Elith & Leathwick, 2009; Barker, Cumming & Darveau, 2014; see Yen, Huettmann & Cooke, 2004 for an application). SAMs are still less commonly used, despite providing valuable information for conservation and management. However, increasing attention has been paid to these problems in recent years (e.g., Yen, Huettmann & Cooke, 2004; Martín et al., 2012; Howard et al., 2015; Ashcroft et al., 2017; Fox et al., 2017).

In the past, spatial conservation decisions and plans were usually just based on SDMs (e.g., Suárez-Seoane et al., 2008; Gray et al., 2009; Adams et al., 2016; Mi, Huettmann & Guo, 2016). However, despite statements by Newton (2008), many scholars found that species occurrence and abundance distribution did not to display similar patterns (Yen, Huettmann & Cooke, 2004; Karlson, Connolly & Hughes, 2011; Yin & He, 2014; Johnston et al., 2015). The difference may represent a mixture of effects and may reflect the differences between the underlying biological processes of abundance and occurrence (Johnston et al., 2015). Therefore, conservation decisions only based on SDMs predictions are insufficient and may even be misleading; the same applies for SAMs. In the future, one time-critical challenge and associated progress will be centered how to combine the useful information that SDMs and SAMs each offer for conservation.

In this study, we evaluated a case study using the endangered Great bustard (Otis tarda dybowskii), which winters in Cangzhou in the North China Plain near Bohai Bay. This area is one of the most important wintering grounds for this species (approximately 300 individuals, c.13.6∼20.0% of China’s total wintering population (Goroshko, 2010; Meng, 2010)). Using the Great Bustard as a case study would contribute to our conservation knowledge about the habitat use of this threatened species and enable us to design better conservation policies. By studying not only the spatial occurrence and the abundance patterns but also combining these two model types together as a role model, predictive modeling and its inferences would potentially have wider conservation implications. Our overall objective of this research was to: (1) assess and develop models to accurately predict the patterns of bustard occurrence and abundance; (2) identify the environmental variables that influence the occurrence and abundance of this species; (3) combine occurrence and abundance models as a new contribution to conservation decisions; and (4) investigate the overall relationship among predicted occurrence, predicted abundance and observed abundance. Well-tested and suitable methods used in this research could be useful for the conservation of the Great Bustard, and other rare species; additionally, this research could generally improve biodiversity through the application of SDMs and SAMs.

Materials and Methods

Study area

This study was conducted at the wintering ground of the endangered Great Bustards in Cangzhou, southeast of the Heibei Province in the wider Bohai Bay (Fig. 1). It is located at 38°12′57″–38°36′51″ latitude and at 116°50′48″–117°24′03″ longitude in the warm temperate, semi-humid monsoon climate zone, which features the slightly marine climatic characteristic of the Bohai Sea region. The topographical and climate conditions vary little in the study area (altitude varies by 13 m, temperature by 0.4 °C, and precipitation is the same). The total study area is 2,191.4 km², consisting of farmland (1,675.1 km²; 76.4%), residential area (330.5 km², 15.1%), open water (23.5 km²; 1.1%) and other unspecified land uses (e.g., home lots, sheds).

Most of the farms in this region produce cereal, which is grown in a 2-year rotation system. In the first year, winter cereal is cultivated from early September to the end of April of the following year. Then, corn is cultivated between the end of April to early September of the same year. The study area was chosen (Fig. 1) because of its large number and proportion (approximately 300 individuals, c.13.6 ∼20.0% of China’s total wintering Great Bustard population (Goroshko, 2010; Meng, 2010)). This area is the world’s largest wintering ground of the endangered O. t. dybowskii. This area is representative of the typical farmland in the North China Plain. In addition, accurate Great Bustard census data, geographic information system (GIS) data coverages and satellite imagery were readily available.

Model development

We employed a machine learning technique, Random Forest, to model the occurrence as well as the abundance distribution of Great Bustards. Breiman’s (2001) Random Forest implementation in SPM7 by Salford Systems Ltd. (San Diego, CA, USA) is robust to over-fitting and is widely recognized to produce high-quality predictive models (Mi et al., 2017). Hence, Random Forest is increasingly applied to species distribution modelling (Cutler et al., 2007; Drew, Wiersma & Huettmann, 2011; Mi, Huettmann & Guo, 2016 for an application using bustards in China). Though Random Forest performed the best in terms of predicting abundance itself (see Appendix S1), testing the feasibility of other data was essential for maintaining high certainty. Thus, to assess the robustness of the model, we pooled data from 2013 and 2014, and then used 80% of the abundance data as training data and the remaining 20% as testing data. When we constructed initial abundance models with all eleven environmental predictors, model performance was not good (the R² value was small). This was likely due to the regression settings in the Random Forest algorithm. For a better outcome, we assessed a “stepwise” setting in SPM for all of the abundance data (100%), re-run the models, and obtained better results. As a result, we identified a multivariate set of four environmental predictors (distance to expressway, distance to national road, distance to pool, and MNNDVI), which had the best performance (the biggest R² value). Using these four predictors, we reconstructed the abundance model based on the training data (80%) and validated it with testing data (20%). We found that the regression model performance was acceptable but fair (R² = 0.551) when comparing observed abundance with predicted abundance. Thus, we constructed the final abundance model based on the above four selected variables described above and the entire set of observation data. To obtain an abundance index more closely aligned with the observations, we adjusted the prediction abundance according to the linear regression between prediction abundance and observation. First, we constructed a regression formula based on the known abundance from observation and prediction [observation abundance = A × prediction abundance + B]. Then, once A and B were be known, they were applied to calculate the adjustment abundance = [adjustment abundance = A × prediction abundance + B]. The regression relation between observation and prediction abundance are shown in Fig. 2.

Further, Random Forest was also applied to rank the relative importance of the environmental variables. In SPMv7, partial dependence plots are not directly implemented in Random Forest but can be obtained easily in R or can be mimicked in the TreeNet model (Friedman, 2002) as a Random Forest run. Thus, we used TreeNet with bagging settings to create partial dependence plots for each variable included in the occurrence and abundance models.

Approximately 10,000 pseudo-absence points were chosen by random sampling across the study areas using the freely available Geospatial Modeling Environment (GME) software (http://www.spatialecology.com/gme/) for distribution models. In SPM7, we set balanced class weights, grew each model to 1,000 classification trees for the occurrence model and 1,000 regression trees for the abundance model, and used all other default setting of the software. We extracted the habitat information for presence and pseudo-absence points for Great Bustards from the environmental layers in GME (“isectpntrst” commands), and then we created a model file in SPM7 called a ‘grove’ that contained the algorithm that quantified the patterns of occurrence; this was used to score all pixels in the study area. We also extracted the habitat information from the same environmental layers for abundance points and then generated a ‘grove’ file for abundance to score abundance estimates for each pixel in the study area.

For spatial occurrence and abundance distribution visualization, we applied the SPM7 grove files to a regular lattice of points (pixels; also attributed to the environmental variables) spaced at 30-m intervals across the study area. Model outputs generated relative indices of occurrence (RIO; an index of pixels from 0 to 1 that represent a relative index belonging to the ‘occurrence’ class) and a relative abundance index (prediction abundance) for each point in the regular lattice based on its underlying environmental variables. We also adjusted the predicted abundance based on a linear regression as constructed in the previous model development steps (Fig. 2A). For a better continuous spatial visualization, the RIO and predicted abundance values were smoothed between neighboring points across the extent of the study area using the Inverse Distance Weighting (IDW) tool in ArcGIS 10.1. This yielded spatially continuous predictive distribution and abundance raster maps of for the Great Bustard.

Results

Partial dependence plots

Partial dependence plots could interpret the functional relationships and effects of each variable by representing a variable’s marginal effects on the response (Elith, Leathwick & Hastie, 2008; Johnstone et al., 2010). It helps to find the signal in the data; Fig. 3A indicates that the occurrence preference of bustards for farmland areas was between 0.6 and 7.5 km². Additionally, based on model results, it appeared the bustard preference were as follows: distance to residential areas ranged from 250 to 2,500 m (Fig. 3B), distance to ditches ranged from 100 to 4,500 m (Fig. 3C), and distance to expressways ranged from 6,000 to 19,000 m (Fig. 3D). In contrast, for abundances, more individuals occurred beyond 2,300 m, but were less than 9,500 m away from national roads (Fig. 3E); additionally, bustards were found in a range between 7,000 and 11,000 m away from expressways (Fig. 3F). Moreover, this species stayed away from pools (maintaining distance greater than 1,500 m, Fig. 3G) and preferred areas with more vegetation (mean NDVI during the investigation was larger than 0.13, Fig. 3H). The information for other variables, which were more marginal, can be found in Appendix S2.

Discussion

The occurrence and abundance models of Great Bustard developed here were designed to identify relevant locations for where to prioritize conservation, and to assess the effects of each variable that influenced this species’ occurrence and abundance (Fig. 3). From a multivariate perspective, the area of farmland, distance to residential areas, distance to ditches and to expressways were among the top four most important predictors for bustard occurrence; in contrast, for the abundance model, the variables consisted of another multivariate package that include the distance to national roads, the distance to expressways, the distance to pools and the mean NDVI (Table 2). We found that high RIO habitats had a fragmented distribution throughout the entire study area (Fig. 4A). The abundance model showed that large populations usually occurred in the central and northwestern parts of our study area (Fig. 4B). The center, north and northeast portions of the study area had a high priority protection index (PI) and had a severely fragmented distribution, including these areas should be the priority sites for protection (Fig. 4C). This not only confirms our own records and, with the help of the PI, can now be quantified and modeled further for more effective conservation application.

In our study area, human disturbance was very strong and represented by indicators, such as density of roads and residential areas (Fig. 1). During our study, we also found other threats to this endangered species: these included farmers grazing their sheep; famers sprinkling poison baits in the wheat fields to present sheep from entering; some bird photographers pursued bustards by walking or following birds while on motor vehicles to take photos, which they wanted to show off to others; hunters with dogs chasing hare and ring-necked pheasant during the day and night; some local people hunting bustards; increasing power lines construction in agricultural lands, resulting in bustards sometimes colliding with wires, getting injured, or even dying, especially flying on foggy days or when in a hurry (Janss & Ferrer, 2000); and the interference of firecracker sounds during the Chinese Spring Festival as well as oil rigs and wind farms. Though carrying a high disturbance can result in stress synthesis (e.g., “death by thousand cuts”), a large number of wintering bustards (approximately 300, c. 13.6 ∼20.0% of China’s total wintering population; Goroshko, 2010; Meng, 2010) still wintered in this area. In times of climate change, it can be assumed the population widens (Mi, Huettmann & Guo, 2016). Thus, this is an area of essential importance for bustards in China, regardless of which perspective is taken. Therefore, a feasible conservation plan should be designed, and based on our model’s prediction result, combined with local public customs and financial support as well as a wider buy-in. In our opinion, improving the local pupulation’s education on animal protection, as has been done over the years, would be useful. The same applies to increasing budgets, enforcement and frequency of patrol by the local management and conservation NGOs in the regions with high PI value, and the local community and government should provide more financial support. However, when designing patrol routes designation in the field, the individuals who are monitoring should avoid getting too close to bustards, so as not to disturb and stress the regular wintering activities of the bustards. For the benefit of this species and its habitats, we suggest not converting crop farmland into nursery farmland; we also encourage farmers to harvest their crops with a machine, which is a more beneficial harvesting method for bustards based on our previous research results (Mi, Huettmann & Guo, 2014). We also highly recommend, if possible, to bury power lines into the ground and to collect hunting guns from the local public.

In this study, occurrence and abundance did not display identical spatial distribution patterns, a result which has been reported in some previous studies (Conlisk, Conlisk & Harte, 2007; Karlson, Connolly & Hughes, 2011; Yin & He, 2014; Johnston et al., 2015). There is actually no reason to assume a presence site represents only one individual animal, nor should a linear relationship between RIO and abundance be assumed. Technically speaking, ‘presence’ can mean one to infinite animals are present, and details depend on the actual pixel setup and how it fits into the obtained model. Therefore, while the relationship is not automatically clear, this could be due to several reasons and depending on specific habitat details: Firstly, the environmental variables that contributed to occurrence and abundance were different, as indicated in Table 2. Secondly, the predictors of preference for bustard occurrence and abundance models were different. For instance, bustards occurred in areas with a distance to expressways ranging from 6,000 to 19,000 m (Fig. 3D), while most populations occurred between 7,000 and 11,000 m from expressways in terms of abundance (Fig. 3F) (see more details in Fig. 3 and Appendix S2). Thirdly, they differed in their spatial distribution for occurrence and abundance (Figs. 4A and 4B). Based on the analysis of overlaying the observation sites with the RIO and observation abundance (Figs. 5A and 5B), the estimated relative index of occurrence (RIO) was not consistently related with the relative index of abundance (Fig. 5A). All locations of observed abundance had high RIO (Fig. 5A), and the relationships between occurrence and abundance estimates were nonlinear (Fig. 5B). These differences may represent a mixture of effects that reflect differences between the underlying biological processes that give rise to specific abundance and occurrence at a specific pixel, as well as limitations imposed by the data and methodology used to estimate these patterns (Johnston et al., 2015; see Buckland et al., 2016 for Distance Sampling and detection problems). In addition, how to interpret the inconsistency between these two indices of plant prediction is a problem waiting to be further resolved e.g., between crop occurrence index (equal to habitat suitability index) and crop abundance (e.g., production).

Treating all presence reports as equal in species distribution models (SDMs:occurrence model, habitat niche model) regardless of the abundance of individuals that the habitat supports could provide us with the information on the loss of habitat suitability (Howard et al., 2014). Applying models based on abundance data, even at a relatively coarse scale, can help to predict spatial patterns of occurrence that are modelled with even greater refinement (Howard et al., 2014). Conservation decision-making should use as much knowledge and information as possible to optimize the benefits of conservation actions (Sutherland et al., 2004; Segan et al., 2011). The use of species distribution models (SDMs) of occurrence has been an important tool in optimizing the selection of protected areas (Franklin, 2013; Guisan et al., 2013; Mi, Huettmann & Guo, 2016; Han et al., 2017) based on the ecological niche space (Drew, Wiersma & Huettmann, 2011). However, relative abundance is often perceived to be a more relevant metric because it can quantify animals within a pixel, and thus, populations (Johnston et al., 2015). Modeling abundance requires methods that can handle large numbers of zero counts as well as the rare, but important, high counts (Welsh et al., 1996), even without a solid research design, according to frequentist statistics. However, Yen, Huettmann & Cooke (2004), Magness, Huettmann & Morton (2008) and Fox et al. (2017) have already shown how machine learning can change this perspective and provide very powerful solutions.

High counts and their locations are particularly important because the pixels with the highest densities of animals are potentially of the greatest interest for conservation planning (Johnston et al., 2015). In our study, we found that the regressions in Random Forest performed poorly at sites with low and high counts (Fig. 2B), although it showed a highly linear relationship between observed and predicted abundance (R² = 0.844; Fig. 2A). Therefore, we argue that the regression method in the Random Forest algorithm should optimize low- and high-count predictions. We recommend to classifying abundances into bins (e.g., high, medium, and low with associated abundance estimates) because Random Forest is exceptionally strong for classification problems. However, for now, this remains an open field of research, but we find our progress remains substantial.

Abundance data could also provide valuable baselines against which to assess future changes (Cumming, 2007) (e.g., climate change, land use change). Such changes in abundance will be much more rapidly apparent, and hence, more rapidly detected than changes in presence-absence patterns across large ranges (Gregory et al., 2005). However, only a few spatial distribution modelers derived models based on the collection of abundance data (e.g., Yen, Huettmann & Cooke, 2004; Fox et al., 2017). This may be because collection of abundance data is more cost or resource demanding than collecting presence—absence data especially for highly mobile animals. Such data are sophisticated in structure and research design, and still they are rarely shared (see in http://GBIF.org). Therefore, we recommend that abundance data could be collected (and can also easily to be turned into presence-absence data, too), even at only relatively coarse numerical scales, because the benefits are considerable (as stated by Howard et al. 2014). One thing that should be mentioned is that plenty of abundance data and models did not perform well and, abundance was extremely difficult to predict (Oppel et al., 2012). Finding the underlying causes that influence abundance model accuracy and constructing more accurate models would be extremely important and useful in future applications toward individual-based policy applications.

For a spatial priority protection of mobile species, one should note that high numbers of individuals are not always present in the same habitats and pixels; instead, low numbers of individuals may occur in one place many times. In addition, this may have implications on spatial priority protection for mobile species. Previous studies have used analytical approaches to address some of these challenges (e.g., Nichols, Thomas & Conn, 2009; Kery & Andrew, 2010; Oppel et al., 2012; Jiguet et al., 2013). However, no general modeling framework has been proposed for dealing with all these analytical challenges simultaneously. This is exactly where our PI offers progress. We also thought the situation of mobile species selecting habitats could be divided into five scenarios: higher numbers and multi-frequency, higher numbers and lower frequency, low numbers and multi-frequency, low numbers and low frequency, and none. When a conservation plan is designed for a species, one should consider not only occurrence index and frequency but also abundance. Here, we proposed the priority protection index (PI; Eq. (1) and Fig. 4), based on the distribution of occurrence and abundance patterns, as a helpful tool for more quickly designs a priority protection plan compared to indices, and it is only based on distribution of occurrence or abundance.

To date, quantitative estimates of population size during global and local changes have actually proven to be difficult to forecast (O’Grady et al., 2004). This is a major hindrance for effective management, as population size and trends are considered among the best correlates of extinction risk (O’Grady et al., 2004). Such measures are commonly used in determining the conservation status of a species (e.g., International Union for Conservation of Nature (IUCN)). We argue that habitat loss remains the one and only powerful metric that can be obtained quickly on a landscape-scale in the absence of proper trend and abundance data (e.g., Drew, Wiersma & Huettmann, 2011). The relationship between predicted environmental suitability and abundance—as presented here—may provide us with a possible method for predicting population size and its changes associated with distributional changes; additionally, it may be particularly appropriate for non-mobile species (e.g., plants, fungi). However, this method is not particularly suitable for mobile species, especially for highly mobile species, such as many birds, bats, and flying insects. They may move over a large landscape within a single day, and abundance and the environment conditions can vary seasonally and spatially. When computing population size or population density using abundance, the primary task will be determining the appropriate unit area for investigation and conservation management.

This study is the first that has combined model-predicted occurrence (representing species distribution models) and abundance indices (representing species abundance models) to produce a priority protection index (PI), which may contribute to spatial conservation and management decisions worldwide. We strongly encourage other researchers to test, apply and update the priority protection index (PI) to explore the generality of these findings further.

Supplemental Information