Prioritizing Green Spaces for Biodiversity Conservation in Beijing Based on Habitat Network Connectivity

: Rapid urbanization results in changes in land use, biogeochemical cycles, climate, hydrosystems, and biodiversity. Policy-makers have formulated ecological protection measures to facilitate sustainable development. However, traditional conservation planning mainly focuses on protecting speciﬁc green spaces, with limited consideration of the connectivity among green spaces from a habitat network perspective. Using citizen science data and occupancy modelling, we predicted habitat suitability, built habitat networks and identiﬁed key habitat patches based on their contribution to the functional connectivity of the habitat network for three focal water, forest, and open-habitat bird species. Based on the habitat requirement, small waterbodies and intermediate forest and open-habitat cover facilitate preserving water, forest and open-habitat birds. In regards to the network analysis, we found that key habitat patches with a high conservation priority were generally characterized by a relatively large patch size and/or located at critical positions in the habitat network (at central positions in the habitat network, or near large patches). We suggest that key habitat patches in restricted built-up areas are converted to protected areas or are kept as cropland under future urban planning. We emphasize the usefulness of the focal species concept in urban biodiversity conservation. Our study offers conservation recommendations from a habitat network perspective for urban planners to safeguard urban biodiversity and ecosystem health.


Introduction
The world has experienced widespread urbanization, especially since 1950, with around 55% of the current human population living in urban areas [1]. In China, the urbanization rate was 18% in 1978, and reached 56% by 2015, due to the program of economic reform starting in 1978 [2]. Rapid urbanization results in changes in land use, biogeochemical cycles, climate, hydrosystems, and biodiversity, from local to regional scales [3]. Furthermore, the conversion of urban green spaces and suburban land into built-up land results in a fragmented landscape and isolated habitat patches for flora and fauna [4]. Such changes negatively affect species' survival, leading to both functional

Bird Data and Focal Species Selection
We downloaded bird observation data from the online gateway China Bird Report Center [43] , which as metadata includes the species name, observation time, observation site (place name), species count and the username of the observer. We used relatively recent bird observation records for the 2008-2017 decade. We checked the location of each site based on 2015 Google Earth images and obtained the longitude and latitude of the site center following Li et al. [44]. In total, we obtained observation records for 545 bird species at 150 observation sites, including 81 resident, 179 migratory, 91 breeding, 42 wintering and 152 vagrant bird species [45,46]. We focus on resident birds (14 generalists were excluded) with observation records available throughout the year. We then categorized them into water, forest, and open-habitat bird communities, according to their main habitat requirement [45,46]. Water, forest, and open-habitat birds refer to groups of birds that use water, forest, or open areas (such as farmland and grassland), as their main habitat.

Study Area
Beijing is the capital of China, located in the northern part of the North China Plain (39.4 • -41.6 • N, 115.7 • -117.4 • E). Beijing municipality covers 16,808 km 2 and is characterized by distinct landforms, with mountains located in the northwest and plains located in the southeast (Figure 1a). Mountainous area occupies 61.3% of the municipality and is mainly covered with shrubs and trees. The plain covers the remaining 38.7% and is the main area of urban sprawl [42]. The diverse habitat types produced by the combination of landforms support a rich flora and fauna.

Bird Data and Focal Species Selection
We downloaded bird observation data from the online gateway China Bird Report Center [43], which as metadata includes the species name, observation time, observation site (place name), species count and the username of the observer. We used relatively recent bird observation records for the 2008-2017 decade. We checked the location of each site based on 2015 Google Earth images and obtained the longitude and latitude of the site center following Li et al. [44].
In total, we obtained observation records for 545 bird species at 150 observation sites, including 81 resident, 179 migratory, 91 breeding, 42 wintering and 152 vagrant bird species [45,46]. We focus on resident birds (14 generalists were excluded) with observation records available throughout the year. We then categorized them into water, forest, and open-habitat bird communities, according to their main habitat requirement [45,46]. Water, forest, and open-habitat birds refer to groups of birds that use water, forest, or open areas (such as farmland and grassland), as their main habitat. The focal species concept is defined as a species or group of species for which their spatial and functional requirements encapsulate those of many other species [47]. For each of the three bird communities we chose one focal species. We used three criteria for selection: the habitat "barycenter" B s [48], the observation frequency (OF), and the variance of B s confidence intervals (VBCI). Habitat "barycenter" B s quantifies the habitat requirement (here referring to the habitat coverage) of a species. VBCI quantifies the variance of the habitat requirement (difference between the upper and lower limit of 95% confidence intervals of habitat coverage at each site). Species with a relatively narrow variance is more selective. Observation frequency (OF) was calculated as the number of sites in which a species was observed, divided by the total number of observation sites. It indicates the sample size of the presence data for each species.
Habitat "barycenter" B s is calculated by Equation (1): where B s is the habitat barycenter of species s, n i is the number of individuals at the ith observation site, x i is the value of environmental variables (water, forest, or open-habitat coverage) at the ith observation site and N is the total number of individuals recorded at all sites.
Since the median area of the observation sites is 1.2 km 2 , which equals the size of a 600 m-radius buffer, we used a 600 m-radius buffer from the center location of the observation site to calculate the habitat composition for each observation site, B s and VBCI. All analyses were carried out using the raster package version 2.8-19 (University of California Davis, Davis, CA, USA) [49] in the R statistical software version 3.5.0 (R Foundation for Statistical Computing, Vienna, Austria) [50].
Species with both high B s (>0.3) [48] and OF (>0.2, equals to being observed at more than 30 sites) values were pre-selected as candidate focal species (Tables A1-A3, Appendix A). We then calculated the weighted overall score (B s × 0.3 + OF × 0.3 + (1 − VBCI) × 0.4) for these preselected species and rank them accordingly. Species with the highest overall score were selected as the focal species, namely mallard (Anas platyrhynchos), grey-headed woodpecker (Picus canus), and Eurasian collared dove (Streptopelia decaocto), as the focal species for water, forest, and open-habitat bird communities.

Environmental Variables
Environmental variables used to predict the habitat suitability for focal water, forest, and open-habitat bird species are summarized in Table 1. Six site-level environmental variables that are ecologically related to habitat suitability were used to predict the occupancy probability. Three observation-level environmental variables were used to quantify observation bias. Land-use composition determines the availability of food and roosting areas [21,57]. We used the best available land cover map, produced in 2015 and with an overall classification accuracy of 88% [51]. Land cover types include water, built-up land, bareland, tree/shrub (forest), grass and cropland. Bareland refers to areas that vegetation is hardly observable but dominated by exposed soil, sand, gravel, and rock backgrounds. We then calculated the percentage of land cover type within a buffer of 600 m around the observation location to extract habitat composition using FRAGSTATS software Version 4 (University of Massachusetts, Amherst, USA) [58].
Additionally, the distance to the nearest water body indicates the accessibility to water resources, another important factor for bird habitat selection [59,60]. Distance to road positively affects bird species richness [61] and traffic volume on roads positively correlated with woodland bird mortality [62]. We then adopted distance to the nearest major road and nightlight index as proxies of human disturbance [61,63,64]. Considering multi-collinearity problems among variables, we kept one variable if two were correlated (a Pearson correlation coefficient >0.5 [65]) and retained the one that was more relevant to bird presence based on expert knowledge. As time of year (date), temperature and precipitation during observations influence detection probability [31], we extracted monthly average temperature and monthly total precipitation for each observation site for the period 2008-2017. All the calculations were conducted in ArcGIS software Version 10.3.1 (ESRI Inc., Redlands, CA, USA) [66].

Habitat Suitability Modelling
Citizen science data are often collected without a standardized field protocol and explicit sampling design [30]. To cope with potential sampling bias, we used occupancy modelling, which takes imperfect detection of species into account and separates detection processes and species' occupancy in a hierarchical way [30,31]. If a species was recorded in a specific month, the detection record was assigned a value of 1 for that month and if not it was assigned a value of 0. In addition, a site not visited at a specific month was considered as a "missing observation". The probability for the observed detection history (0/1) of a species at a specific site P observed is the combination of the occupancy probability ψ and the detection probability (p, Equation (2)): where ψ and p are functions of site-level and observation-level variables, respectively (Equations (3) and (4)).
Logit p ij = B 0 + B 1 y j1 + B 2 y j2 + · · · + B n y jn  [67] in the R statistical software version 3.5.0 (R Foundation for Statistical Computing, Vienna, Austria) [50]. We performed model selection on both occupancy and detection probability submodels simultaneously using the R package MuMIn version 1.42.1 (Institute of Nature Conservation Polish Academy of Sciences, Krakow, Poland) [68] in the R statistical software version 3.5.0 (R Foundation for Statistical Computing, Vienna, Austria) [50]. We used Akaike's Information Criterion adjusted for small sample size (AICc) to rank the models ( AICc < 2) [69]. We examined the goodness of fit of selected models using a chi-square (χ 2 ) test with 1000 bootstraps replicates and excluded models with a p-value larger than 0.05 [70]. Considering some top models might have a lower Akaike weight, we adopted a model-averaging approach [69] by starting with the model with the highest Akaike weight and including the model with the second highest weight until ∑ w m = 0.95. Then the coefficient of each variable in the final model was estimated by the weighted mean of the selected models (according to the Akaike weight of the model) [71]. The response curve of the most important environmental factors for each focal species was generated by plotting the response value (and the 95% confidence intervals) at the range of the variable of interest, based on the prediction of the final occupancy submodel.
To define suitable and unsuitable area based on the predicted occupancy probability (0-1), we used the sample frequency instead of an arbitrary number as a threshold [72]. Sample frequency was calculated as the proportion of species' occurrences across all sites [72]. Hence, areas with an occupancy probability above 0.55, 0.57, and 0.79 were defined as suitable habitat to build the habitat network for water, forest and open-habitat focal bird species. Generating habitat suitability maps and extracting suitable habitat patches were conducted in ArcGIS software [66].

Habitat Network Connectivity Analysis
The suitable habitat patches, as stepping-stone nodes, are crucial connections to maintain habitat network connectivity. We employed the probability of connectivity (PC), which is considered a robust metric to evaluate the functional connectivity of a habitat network [15,20]. If there is a set of habitat patches and connections among them, PC is defined as the probability that two animals are randomly placed within the landscape and fall into habitat areas that are mutually reachable [15], and is calculated with Equation (5): where n is the total number of habitat patches, a i and a j represent the area of suitable habitat patch i and j, p * ij is the maximum product probability of all the possible paths between patches i and j, and A L stands for the total study area.
We quantified the importance of patches by removing each single patch from the network to determine the change in network connectivity [15,20]. The importance of each patch is assessed by the absolute difference in the probability connectivity (dPC) before and after the removal, and is calculated as a percentage from Equation (6): where PC k is the overall connectivity value when all existing patches are present in the landscape and PC remove,k is the connectivity value after the removal of a specific patch k. A patch with a greater dPC has a higher contribution to the network connectivity. Different bird species vary in their dispersal abilities. Hence, we performed connectivity analysis over a range of dispersal distances from 1-50 km (1, 3, 5, 10, 20, 30, 40, 50 km), and adopted the mean of the eight dPC values to make the results applicable to birds with different dispersal distances. The habitat network connectivity analyses were performed using the Conefor Inputs extension in ArcGIS software [66] and Conefor Sensinode software Version 2.2 (University of Lleida, Lleida, Spain) [73]. Based on their contribution to habitat network connectivity (dPC), we classified habitat patches into five levels of conservation priority (very high, high, medium, low, and very low) using the natural break (Jenks) method [74] in ArcGIS (grouping similar dPC values into one class and maximizing differences between classes). Patches assigned with a high or very high level of conservation priority were defined as key habitat patches. We then identified whether key patches are protected or not, or whether they fall into the restricted build-up areas, based on the boundary of the ecological control area outlined in the Beijing Urban Master Plan (Figure 1b).

Habitat Suitability for Focal Species
The final models for the three focal species are shown in Table 2. Mean detection probability of the mallard (water bird) was the highest (0.56), followed by the grey-headed woodpecker (forest bird) (0.40), and then the Eurasian collared dove (open habitat bird) (0.20). Table 2. Final models for estimating the occupancy (ψ) and detection (p) probability for focal water bird mallard, forest bird grey-headed woodpecker and open-habitat bird Eurasian collared dove in Beijing.
For mallard, the percentages of water and grass were positively correlated with occupancy probability ( Table 3). The most important variable was water and mallard only required relatively low coverage of water (Figure 2a). Suitable habitats for water birds mainly appeared around reservoirs, wetlands, lakes, and rivers (Figure 3a), such as Miyun reservoir in Northeast Beijing, Guanting reservoir and Wild Duck Lake national wetland park in the northwest, as well as wetlands and green spaces along the Chaobai and Wenyu rivers in urban areas (Figure 3a). For grey-headed woodpecker, the percentage of forest positively influenced occupancy probability and was also the most important variable ( Table 3). The occurrence probability of this species sharply increased when the forest cover reached 30% (Figure 2b). Highly suitable areas were concentrated in the northern and southwestern suburban mountain forest and forest in urban green spaces (Figure 3b). For Eurasian collared dove, the percentage of open habitat (grass and cropland) positively influence occupancy probability ( Table 3). The probability of presence increased considerably when the cover of grass and cropland reached 25% (Figure 2c). Suitable habitats for open habitat birds were mainly positioned in Eastern and Southern Beijing in the Fangshan, Daxing, Tongzhou, Shunyi, Pinggu, and Yanqing districts (Figure 3c).    Open habitat refers to grassland and cropland. The grey areas indicate the 95% confidence interval.

The Habitat Network Connectivity and Conservation Priority
For the focal water bird species, key habitat patches were mainly distributed in east and central Beijing (the Chaobai river in Shunyi district, the Liangshui and Xiaotaihou rivers in Tongzhou district, the Wenyu and Ba rivers in Chaoyang district, Miyun reservoir in Miyun district, the North Canal and Yunchaojian rivers in Tongzhou district, the Ju and Ru rivers in Pinggu district, the Sha and Wenyu rivers in Changping district, the Qing river in Haidian district and Olympic Park in Chaoayng district) (Figure 4a   Most of the key habitat patches for forest birds (92%) were located in ecological control areas, whereas a small percentage fell into restricted built-up areas ( Figure 5). For water and open-area birds, only 37% and 52% of key habitat patches were in ecological control areas ( Figure 5). Key habitat patches for water birds in restricted built-up areas, with a total surface area of 321.3 km 2 , were mainly distributed in Northeast Beijing (along the Chaobai, Ru, and Ju rivers), Central Beijing (surroundings of Shahe reservoir and the Wenyu river), and Southwest Beijing (along the Liangshui river in Yizhuang). For forest birds, only 272.9 km 2 of key habitat patches were located in restricted built-up areas, and mainly located at the foot of the mountains and in the Xiaotangshan modern agricultural park in Changping district. A total of 1036.5 km 2 of key habitat patches for open-habitat birds occurred in restricted built-up areas, and were mainly distributed in the Daxing, Tongzhou, and Shunyi districts.

Discussion
Our study proposes an approach in prioritizing urban green space protection from a focal species habitat network perspective. Using occupancy modelling and citizen science bird observation data, we predicted habitat suitability for three focal water, forest, and open-habitat bird species,

Discussion
Our study proposes an approach in prioritizing urban green space protection from a focal species habitat network perspective. Using occupancy modelling and citizen science bird observation data, we predicted habitat suitability for three focal water, forest, and open-habitat bird species, presumed to encapsulate the needs of wider species communities. Since different communities show distinct habitat preference, our approach could more accurately quantify habitat patches for targeting communities, compared to simply using green space to delineate the habitat in general. The identified core habitat patches in the habitat networks could guide urban green space management and provide the scientific basis for urban planners to target areas for protection priority in the future development plan. By using focal species of three bird communities as bio-indicators, our study provides evidence-based measures to help maintain urban habitat connectivity and preserve biodiversity.
The quantified habitat requirements, suitable habitat maps and identified core habitat patches for the focal water, forest and open-habitat bird species could guide green space management and urban planning. The focal water bird mallard requires only a relatively low surface area of water, suggesting maintaining or restoring small waterbodies in urban green areas would be an efficient way of maintaining water birds. An alternative explanation would be that most waterbodies in urban areas are relatively small and water birds have adapted to this kind of landscape. On the other hand, frost birds and open-habitat birds need a moderate level of forest and open-habitat cover to survive. This implies that the composition of forest and open-habitat should be taken into account when designing city parks, restoring habitats or setting up conservation plans.
We found that large habitat patches are important in facilitating habitat network connectivity. For example, patch 'B', 'C', and 'G' (Figure 4) are the largest patches in the habitat network of the focal species and are also some of the most important patches in facilitating habitat connectivity. Patch 'B' is located in an area composed of rivers (e.g., the Liangshui and Xiaotaihou rivers), wetlands, parks and grassland. Patch 'C' is covered by continuous forest and 'G' is grass/farmland. Larger patch size contributes to a greater intra-patch and inter-patch connectivity in a habitat network [16,19]. They not only offer a larger area of habitat, but also increase the amount of area-weighted dispersal flux [19] and, hence, should to be protected.
In addition to the size of patches, we found that habitat patches positioned at a critical spatial position in the habitat network are also important for maintaining habitat connectivity. By acting as stepping-stones, patches located in the center of the network and/or near large patches, tend to have a higher contribution to the habitat connectivity. For example, for the representative water bird species, patch 'A' (Figure 4a) along the Chaobai River connects surrounding habitat patches. For the representative forest bird species, patch 'D' (Figure 4b) along the Wenyu River, acted as a key stepping stone connecting the urban and suburban areas. For the focal open-habitat bird species, patch 'E' (Figure 4c) connects the northern and southern key patches. Moreover, patches near large patches (e.g., patch 'F' in Figure 4c) not only facilitates species dispersal, but also acts as a barrier against human disturbance [16]. Loss of such key stepping-stone patches might rapidly reduce the overall connectivity and increase the chance of network collapse [75,76]. It is, therefore, important to evaluate sites for conservation priority from a network perspective, instead of only at the site level. Hence, such stepping-stone patches should be prioritized in urban biodiversity conservation [16,37].
We suggest conservation planning takes into account both habitat suitability and habitat network connectivity. As the designation for restricted built-up areas in Beijing is not yet fixed, we suggest that the identified key habitat patches within restricted built-up areas for water and forest birds should be converted into ecological control areas and those for open-habitat communities should be preserved as cropland under future urban planning (Table A4, Appendix A).
We selected focal water, forest and open-habitat birds with distinct habitat preferences. The spatial and functional requirement of these focal species could help define the environmental limits and help protect other species [33]. A rigid selection procedure is critical, given that focal species need to reflect the overall biodiversity of the territory under investigation. Besides the commonly used criteria habitat barycenter [33,48], we also considered the variance of the habitat barycenter. Smaller variance indicates that the species is more selective toward one specific habitat type and uses alternative habitats less. In consequence, by protecting the habitat of selected focal species, other taxonomic groups could also be protected, including those with less observation data available. We underline the usefulness of applying the focal species concept in urban biodiversity conservation.
Our study highlights the combination of occupancy modelling of focal species with network connectivity analysis to identify key habitat patches for conservation prioritization. Species distribution modelling (SDM) has been frequently used to quantify environment-species relationships and generate suitability habitat maps [77] to guide urban wildlife protection and biodiversity conservation (e.g., [32,78]). The habitat network connectivity has hardly been investigated after suitable habitat has been delineated using SDM (except for [22,48,79]). In real world conservation, not all suitable habitat can be protected, especially under ongoing urbanization. The quantification of the core habitat patches from a habitat network connectivity perspective offers an efficient way to target the most important areas for protection priority. In addition, we built the occupancy model based on 10-year bird observation data and the 2015 land cover map, mainly due to the limited number of bird records. We recommend promoting the collection of citizen science data, as more observations and observation sites would improve the prediction of suitable habitat patches. Nevertheless, our approach of coupling focal specie occupancy modelling with network analysis to guide conservation efforts could be applied to other cities, where fine resolution land cover map and bird observation data are available.

Conclusions
Since habitat fragmentation causes biodiversity loss, it is necessary to protect urban biodiversity based on habitat suitability and the functional connectivity of the habitat network. We select three focal species presuming their habitat requirements encapsulate many other species and underline their usefulness in urban biodiversity conservation. Based on the habitat requirement, small waterbodies and intermediate forest and open-habitat cover facilitate preserving water, forest, and open-habitat bird communities, respectively. In light of the network analysis, key habitat patches with a high conservation priority are generally characterized by having a large surface area and/or being located at critical positions in the habitat network (i.e., at central positions in the habitat network and/or near large patches). We suggest that urban planners preserve such key patches in restricted built-up areas (Table A4, Appendix A) in future urban development plans. We consider that the evaluation of individual habitat patches in terms of conservation priority, based on the contribution of these patches to habitat network connectivity, is an effective method to optimize decision-making in future urban conservation planning. Our findings provide guidelines for the conservation of local biodiversity and sustainable urban development.