GIS-Based Assessment of Habitat Networks for Conservation Planning in Kas-Kekova Protected Area (Turkey)

The determination of protected area (PA) boundaries and the level of restrictions is very important for sustainable conservation, and such decisions must involve biodiversity data and assessment. In a PA, the consensus of the government and the locals is crucial to ensure protection–use balance. The PA restrictions constrain legal human activities, and the boundary determination and the restrictions should be based on various scientific analyses to achieve consensus. In this study, a GIS-based approach is proposed to utilize the biodiversity data for efficient conservation and land use planning in Kas-Kekova PA, which is among the most important PAs in Turkey. Spatial analysis methods, i.e., kernel density estimation, natural breaks classification and integrated density index, were performed for the assessment of the habitat networks using georeferenced biodiversity datasets, and the results were evaluated with respect to the actual land use data and the land ownership pattern. The developed spatial analysis approach is efficient to produce the conservation base maps required for regional land use planning, for defining sustainable conservation strategies, and to provide a widely accepted base for land use planning and biodiversity monitoring in the PA; although careful investigations and expert opinions are still required for data deficient areas.


Introduction
In the Earth Summit on Development carried out in Rio de Janeiro in 1992, it was accepted that the loss of biodiversity is an important issue that can be prevented by coordinated global efforts. The Convention on Biological Diversity (CBD) was adopted for "sustainable development" in the same year [1]. The Convention mainly focused on the conservation of biological diversity. Despite all efforts, the loss of biodiversity continues, and both climate change and environmental pollution play significant roles in this process [2][3][4]. Therefore, integrated global efforts should be increased to protect natural habitats and to preserve them for future generations.
Many species are in danger of extinction on Earth. The importance of nature protection is emphasized by many national and international organizations. As a leading global environmental authority, the United Nations Environment Program (UNEP) includes the conservation and sustainable use of biodiversity in the global environmental agenda. One of the global biodiversity conservation goals, Aichi Target 11, defined in the CBD strategic plan 2011-2020, has the aim that "By 2020, at least 17% of terrestrial and inland water areas and 10% of coastal and marine areas, especially areas of particular importance for biodiversity and ecosystem services, are conserved through effectively and equitably managed, ecologically representative and well-connected systems of protected areas and other effective area-based conservation measures, and integrated into the wider landscapes and seascapes." [5]. According to [6], there were a total of 202,467 terrestrial and inland water protected areas (PAs) recorded in the World Database on Protected Areas (WDPA) by 2016, covering 14.7% (~20 million km 2 ) of the world's extent of these ecosystems. Although the value falls short of Aichi Target 11, the complexity of the problem and the amount of the geospatial data involved in the processes suggest efficient utilization of geographical information systems (GIS) and related technologies.
Designation of a conservation area as PA by governments does not stop the loss of biodiversity by itself [5]. There are many issues that should be organized in a PA to implement a successful plan [7,8]. In addition, strict regulations based on ecological networks must be determined [9][10][11]. Essentially, effective protection can only be performed at a local scale. Environmental planning helps to determine the suitable land use by clarifying the biological, physical and social systems of a PA [12]. The first step to understand ecological processes is to identify the ecological patterns. For this purpose, the spatial distribution and densities of flora and fauna elements can be analyzed using spatial statistics [13]. For example, plant species with sparse distributions can be expressed with discrete points [14]. As output of spatial analysis, habitat patterns can be obtained and assessed for defining the regional ecological network. Such a network can be used for improving the habitat connections and proposed as a solution to the land fragmentation problems while ensuring the protection of threatened species and biodiversity [15,16]. As part of an ecological network, ecological corridors and their preservation are also important for sustainable protection. The determination of ecological corridors has important objectives, such as creating a scattered habitat, compensating for habitat losses, and ensuring the protection of endangered species and populations [15].
Several methods were proposed in the literature to identify ecological networks using spatial analysis methods. Opdam et al. [17] proposed a spatial indexing approach to analyze the conservation potential of landscapes for different species. McHugh and Thompson [18] studied an ecological network model for conservation planning using eco-profiles of species. Ferretti and Pomarico proposed an effective tool for spatial planning by integrating GIS with multi-criteria analysis [19]. Gurrutxaga et al. [20] proposed a least-cost path analysis for the determination of ecological corridors based on core habitat zones. Guo and Liu [21] and Santos et al. [22] also proposed an ecological network determination approach based on least-cost analysis at regional scales. Chang et al. [23] proposed a green infrastructure planning approach-based ecological connectivity assessment in Longgang District of Shenzhen in China. Xie et al. [24] proposed an integrated index which evaluates the biodiversity conservation and the disaster protection at the same time. Hong et al. [25] also developed an index system based on ecological sensitivity analysis. The kernel density estimation (KDE) method has been employed very often in the literature to assess the spatial patterns of habitats (e.g., [26][27][28][29]); although its use for evaluating other discrete point pattern analysis applications are also available (e.g., [30,31]). However, the determination of ecological networks is a complex problem and an active research field.
In this study, a GIS-based approach for the determination of habitat networks in Kas-Kekova PA, which is located in the south of Turkey, was investigated. The PA contains diverse endemic flora and threatened fauna species. In addition, there are archeological PAs inside the region. The PA is under threat from biodiversity deterioration or disappearance due to the anthropogenic pressures such as tourism and rapid urbanization. In the area, the use of resources by humans should be managed through sustainable planning and, at the same time, their interest in protecting the biodiversity and the environment should be ensured. An integrated planning approach can assist the development and implementation of conservation principles and support the needs of locals. Therefore, the main goals of the present study are; (a) to develop an easy-to-use approach for producing the base maps to present and involve the biodiversity data for defining land use restrictions in Kas-Kekova PA; (b) to integrate the threatened biodiversity elements (discretely collected flora and fauna data) in land use planning; (c) and to propose additional conservation principles on top of the legal principles defined by the Ministry of Environment and Urbanization (MoEU) based on the spatial analysis results of the biodiversity data, the actual land use and the property ownership pattern. In addition, utilizing a GIS platform to store and present diverse data can help to increase the public and government awareness and to support the decision making.
The article is structured as follows. The main characteristics of Kas-Kekova PA are described in Section 2. The spatial analysis methods and the datasets are provided in Section 3. In Section 4, the results are provided and integrated with the land ownership and land use data, and conservation strategies for land use planning are proposed. The main findings of the study are discussed in the final section.

Main Characteristics of the Kas-Kekova Protected Area
Kas-Kekova Region is located on the coast of Antalya Province, Turkey, between the towns of Kas in the west and Kale in the east (Figure 1). The study area covers a total area of approximately 260 km 2 (~100 km 2 terrestrial and~160 km 2 marine parts). Three different PAs have been declared in the region by two governmental organizations (MoEU and Ministry of Culture and Tourism) based on different laws (Articles 644 and 2863 in [32,33]). The study area was declared as a Highest (First) Degree Natural Protection Site by MoEU in 1989, then as a PA in 1990 under the UNEP and the Barcelona Convention [34]. The PA also contains archaeological ruins, mostly from Lycians, and coastal landscapes that are subject to protection ( Figure 1). principles on top of the legal principles defined by the Ministry of Environment and Urbanization (MoEU) based on the spatial analysis results of the biodiversity data, the actual land use and the property ownership pattern. In addition, utilizing a GIS platform to store and present diverse data can help to increase the public and government awareness and to support the decision making. The article is structured as follows. The main characteristics of Kas-Kekova PA are described in Section 2. The spatial analysis methods and the datasets are provided in Section 3. In Section 4, the results are provided and integrated with the land ownership and land use data, and conservation strategies for land use planning are proposed. The main findings of the study are discussed in the final section.

Main Characteristics of the Kas-Kekova Protected Area
Kas-Kekova Region is located on the coast of Antalya Province, Turkey, between the towns of Kas in the west and Kale in the east (Figure 1). The study area covers a total area of approximately 260 km 2 (~100 km 2 terrestrial and ~160 km 2 marine parts). Three different PAs have been declared in the region by two governmental organizations (MoEU and Ministry of Culture and Tourism) based on different laws (Articles 644 and 2863 in [32][33]). The study area was declared as a Highest (First) Degree Natural Protection Site by MoEU in 1989, then as a PA in 1990 under the UNEP and the Barcelona Convention [34]. The PA also contains archaeological ruins, mostly from Lycians, and coastal landscapes that are subject to protection ( The Kas-Kekova PA Biological Diversity Project was carried out in the study area between 2008 and 2010 by the MoEU [36]. Within the project, European Union Nature Information System (EUNIS) habitat classes, plant species, mammalian, amphibian, reptile and bird species in the PA were determined according to the EU Habitats Directive [37] and EU Birds Directive [38]. The extensive field works held by specialists throughout the project (four consecutive spring and fall seasons in two years) ensured the representative data collection for biodiversity. The PA also has high accessibility by foot and, thus, the whole area could be wandered by the ecology experts in the project time frame. In order to cover the whole PA, the region was divided into 35 square grids with 2.5 × 2.5 km grid size and each grid was numbered [36]. During the eight field campaigns between October 2008 and June 2010, the habitat classes were assessed according to EUNIS. In addition, the field work results were supported with the satellite imagery for accurate delineation of vegetation class boundaries. The endemic and threatened flora were identified during the field works and their locations were measured with ca. 1 m positional accuracy using hand-held GPS (global positioning system) devices. Although mammals and reptiles were observed by the experts, their locations have a lower spatial accuracy since the measured locations were mostly the positions of the observer with an approximated correction. Still, a spatial accuracy of better than 10 m can be expected. The birds' locations have bigger measurement ambiguity, since the identifications were performed by the observers and by voice recordings, and the positions mostly represent the locations of the recorders and observes. Thus, for the birds, a spatial accuracy of better than 50 m can be expected.
A total of 51 families, 187 genera and 272 flora species were identified in the terrestrial habitats of the Kas-Kekova PA. Among those, 26 species are endemic to Turkey. The species were also assessed according to the International Union for Conservation of Nature (IUCN) Red List Categories and Criteria [39,40], which is an important indicator of the sustainability of the world's biodiversity. The IUCN categories of threat are listed as; extinct (EX), extinct in the wild (EW), critically endangered (CR), endangered (EN), vulnerable (VU), near threatened (NT), and least concern (LC). The species in the PA fall into the categories of EN, VU, NT, and LC. Those flora and fauna species which were classified as CR, EN, VU and NT (Table 1) were employed in the spatial analysis algorithms in this study, since they are threatened and their habitats are subject to protection. The reptiles and amphibian data were grouped and mentioned as reptiles in the following sections, since few amphibians were identified in a very small area in the PA.
Kas-Kekova PA also has four major plant communities [36,41], which are pine forests (Aetheorhizo Bulbosae-Pinetum Brutiae), maquis (Quercus aucheri-Oleetum europaeae), phrygana (Alysso-Genistetum acanthocladae) and halophyte (Salicornietum ramosissimae). The maquis community is the main plant association in the PA [36]. This community has a rich and stable biodiversity, includes endemic species, and provides the climax phase and the floristic composition in the PA. The phrygana community was developed as a result of the burning of the maize vegetation in the region and is represented in smaller parts of the PA. The halophyte community develops only in the salt marsh and a very small halophyte plant region exists in the PA, which occurred by the tidal movements coming from the sea. The forest community is mainly composed of red pine and was determined as the union of Aetheorhizo bulbosae and Pinetum brutiae although the Pinus brutia is very rare. These communities were integrated into the land use map, which was generated according to the EUNIS classifications. There are also anthropogenic affects in Kas-Kekova and several rural settlements in and around the PA exist ( Figure 2). The human population inside the PA is ca. 2000 (as of year 2013). Since the PA has a stunning coastal landscape, the tourism activities in the summer season are high and double the population. The daily tourism activities are mostly in the bays and may expand to spring seasons as well. The main income sources of the locals are agriculture, tourism and fishing [36]. Since the 1990s, touristic accommodations and roads in the PA have increased rapidly ( Figure 3). Due to improper infrastructure planning, marinas and hotels prompt illegal construction, which causes the uncontrolled use of natural resources [36]. Although the PA has a master land use plan with a scale of 1:25.000 from 1991, this plan was not based on biodiversity data and does not comply with the human needs either, despite the fact that 20% of the cadastral parcels were registered under private ownership.
Since the economic expectations are more important than environmental protection for locals [42], the conservation strategies must be developed to obtain a protection-use balance.

Materials and Methods
In order to integrate the threatened biodiversity elements into the final land use plans, the spatial habitat patterns of these elements (i.e., endemic and threatened flora and fauna species) should be revealed. This task can be achieved with spatial analysis methods although several complexities should be considered, such as; i) the biodiversity is high (i.e., many species types exist) and performing extensive habitat analysis for every species type in the PA is extremely difficult; ii) only a small number of observations for the threatened species is available and the habitat pattern may not be revealed adequately; and iii) there are large variations between the numbers of observations for different species types. Various spatial analysis approaches can be utilized for habitat pattern determination purposes, but special care must be taken when selecting the method and its parameters. For rarely observed and threatened fauna species, a home range estimation method can be used but it may require exhaustive field observations over a long term. If the number of observations is sufficient for demonstrating the species spatial distributions, a statistical density estimation method can be employed by considering the distribution pattern (e.g., random or clustered). In this study, an easy-to-use methodology was developed for the spatial habitat pattern determination to be used as a base for land use planning in the PA, where the different elements were

Materials and Methods
In order to integrate the threatened biodiversity elements into the final land use plans, the spatial habitat patterns of these elements (i.e., endemic and threatened flora and fauna species) should be revealed. This task can be achieved with spatial analysis methods although several complexities should be considered, such as; i) the biodiversity is high (i.e., many species types exist) and performing extensive habitat analysis for every species type in the PA is extremely difficult; ii) only a small number of observations for the threatened species is available and the habitat pattern may not be revealed adequately; and iii) there are large variations between the numbers of observations for different species types. Various spatial analysis approaches can be utilized for habitat pattern determination purposes, but special care must be taken when selecting the method and its parameters. For rarely observed and threatened fauna species, a home range estimation method can be used but it may require exhaustive field observations over a long term. If the number of observations is sufficient for demonstrating the species spatial distributions, a statistical density estimation method can be employed by considering the distribution pattern (e.g., random or clustered). In this study, an easy-to-use methodology was developed for the spatial habitat pattern determination to be used as a base for land use planning in the PA, where the different elements were

Materials and Methods
In order to integrate the threatened biodiversity elements into the final land use plans, the spatial habitat patterns of these elements (i.e., endemic and threatened flora and fauna species) should be revealed. This task can be achieved with spatial analysis methods although several complexities should be considered, such as; (i) the biodiversity is high (i.e., many species types exist) and performing extensive habitat analysis for every species type in the PA is extremely difficult; (ii) only a small number of observations for the threatened species is available and the habitat pattern may not be revealed adequately; and (iii) there are large variations between the numbers of observations for different species types. Various spatial analysis approaches can be utilized for habitat pattern determination purposes, but special care must be taken when selecting the method and its parameters. For rarely observed and threatened fauna species, a home range estimation method can be used but it may require exhaustive field observations over a long term. If the number of observations is sufficient for demonstrating the species spatial distributions, a statistical density estimation method can be employed by considering the distribution pattern (e.g., random or clustered). In this study, an easy-to-use methodology was developed for the spatial habitat pattern determination to be used as a base for land use planning in the PA, where the different elements were grouped under four categories (i.e., mammals, reptiles, birds and endemic flora) and every observation in each group was represented in the final map properly (e.g., for the threatened mammals) (Figure 4). grouped under four categories (i.e., mammals, reptiles, birds and endemic flora) and every observation in each group was represented in the final map properly (e.g., for the threatened mammals) (Figure 4). First, the flora and fauna point data collected in the fieldwork [36] were analyzed to determine the patterns of plant and animal ecology and to produce ecological sensitivity maps [43,44]. After the initial investigations on the biodiversity reports (i.e., endemic species, habitat types, plant community, threatened categories, location etc.), a geodatabase was populated with spatial and attribute data. The flora and fauna points were filtered for endemism and/or threatened category (EN, VU, NT) attributes and divided into four groups; namely, flora, mammals, reptiles (including amphibians), and birds. After the visual investigations on the filtered point locations, overlaid with the EUNIS habitat map and the topography, a spatial methodology was developed. The KDE method was used for defining the point patterns and an optimal bandwidth determination approach was developed according to the spatial statistics of the points. The natural breaks classification was used to obtain the spatial configurations of each group. The integrated density index was used to combine the patterns of all four groups and thus to obtain a habitat network. The final results were analyzed together with the EUNIS land use classes, land ownership data (cadastral data), and the roads. Several conservation strategies for the PA management are proposed based on the integrated analysis results. The spatial analyses were performed using ArcGIS 10.2 software from ESRI, Redlands, CA, U.S.A. [45].

Spatial Datasets
Kas-Kekova PA has fourteen EUNIS habitat classes ( Table 2) defined according to EUNIS Habitat Classification [46]. The class boundaries were delineated manually using satellite images and inspected by field work [36]. The formats, sources, and the accuracies of datasets evaluated in the study are given in Table 3. Figure 5 shows the EUNIS Habitat Map of the PA and the locations of threatened fauna and endemic flora. The region is covered largely by forests (82.4%). In Figure 6, the cadastral parcels and their ownerships are demonstrated. The ownerships are categorized as; i) private ownership; ii) private ownership status obtained after deforestation in the area, iii) parcels owned by the state, and iv) state forest areas that are not subject for private ownership. First, the flora and fauna point data collected in the fieldwork [36] were analyzed to determine the patterns of plant and animal ecology and to produce ecological sensitivity maps [43,44]. After the initial investigations on the biodiversity reports (i.e., endemic species, habitat types, plant community, threatened categories, location etc.), a geodatabase was populated with spatial and attribute data. The flora and fauna points were filtered for endemism and/or threatened category (EN, VU, NT) attributes and divided into four groups; namely, flora, mammals, reptiles (including amphibians), and birds. After the visual investigations on the filtered point locations, overlaid with the EUNIS habitat map and the topography, a spatial methodology was developed. The KDE method was used for defining the point patterns and an optimal bandwidth determination approach was developed according to the spatial statistics of the points. The natural breaks classification was used to obtain the spatial configurations of each group. The integrated density index was used to combine the patterns of all four groups and thus to obtain a habitat network. The final results were analyzed together with the EUNIS land use classes, land ownership data (cadastral data), and the roads. Several conservation strategies for the PA management are proposed based on the integrated analysis results. The spatial analyses were performed using ArcGIS 10.2 software from ESRI, Redlands, CA, U.S.A. [45].

Spatial Datasets
Kas-Kekova PA has fourteen EUNIS habitat classes ( Table 2) defined according to EUNIS Habitat Classification [46]. The class boundaries were delineated manually using satellite images and inspected by field work [36]. The formats, sources, and the accuracies of datasets evaluated in the study are given in Table 3. Figure 5 shows the EUNIS Habitat Map of the PA and the locations of threatened fauna and endemic flora. The region is covered largely by forests (82.4%). In Figure 6, the cadastral parcels and their ownerships are demonstrated. The ownerships are categorized as; (i) private ownership; (ii) private ownership status obtained after deforestation in the area, (iii) parcels owned by the state, and (iv) state forest areas that are not subject for private ownership.       The topographical characteristics of the PA were analyzed in terms of slope, aspect and elevation values and investigated with respect to mammals, reptiles and endemic flora species. The slope values in the PA range between 0-33 degrees (Figures 7 and 8) and the elevation values range between 0-1170 meters. Although topography may form barrier to biodiversity elements, considering the spatial distributions of the mammals (i.e., Capra aegagrus Erxleben) and the reptiles locations identified in the PA, the topography was not found to be an effective factor or barrier (Figures 7 and 8). Therefore, it was not considered in the ecological network assessment. In addition, no clear indications of spatial correlation with the topographical characteristics (slope, aspect and elevation ranges) and the flora data were observed as shown in Figure 9.  The topographical characteristics of the PA were analyzed in terms of slope, aspect and elevation values and investigated with respect to mammals, reptiles and endemic flora species. The slope values in the PA range between 0-33 degrees (Figures 7 and 8) and the elevation values range between 0-1170 m. Although topography may form barrier to biodiversity elements, considering the spatial distributions of the mammals (i.e., Capra aegagrus Erxleben) and the reptiles locations identified in the PA, the topography was not found to be an effective factor or barrier (Figures 7 and 8). Therefore, it was not considered in the ecological network assessment. In addition, no clear indications of spatial correlation with the topographical characteristics (slope, aspect and elevation ranges) and the flora data were observed as shown in Figure 9. The topographical characteristics of the PA were analyzed in terms of slope, aspect and elevation values and investigated with respect to mammals, reptiles and endemic flora species. The slope values in the PA range between 0-33 degrees (Figures 7 and 8) and the elevation values range between 0-1170 meters. Although topography may form barrier to biodiversity elements, considering the spatial distributions of the mammals (i.e., Capra aegagrus Erxleben) and the reptiles locations identified in the PA, the topography was not found to be an effective factor or barrier (Figures 7 and 8). Therefore, it was not considered in the ecological network assessment. In addition, no clear indications of spatial correlation with the topographical characteristics (slope, aspect and elevation ranges) and the flora data were observed as shown in Figure 9.    (Table 1) observed during the field works on the altitude range (a), slope (b) and aspect (c) maps.

Spatial Analysis Methods
The KDE method has often been used in the literature for the density estimation of point patterns (e.g., [26][27][28][29][30][31]). It is a statistical method developed by Silverman [47] and acts as a smoothing filter for discrete points (Equation 1). The KDE method is based on the value of i at each sampling point (x, y) in a kernel. The density contribution of each sampling point was calculated using the kernel functions. The density value changes depending on the search radius range (bandwidth) and grid   (Table 1) observed during the field works on the altitude range (a), slope (b) and aspect (c) maps.

Spatial Analysis Methods
The KDE method has often been used in the literature for the density estimation of point patterns (e.g., [26][27][28][29][30][31]). It is a statistical method developed by Silverman [47] and acts as a smoothing filter for discrete points (Equation 1). The KDE method is based on the value of i at each sampling point (x, y) in a kernel. The density contribution of each sampling point was calculated using the kernel functions. The density value changes depending on the search radius range (bandwidth) and grid  (Table 1) observed during the field works on the altitude range (a), slope (b) and aspect (c) maps.

Spatial Analysis Methods
The KDE method has often been used in the literature for the density estimation of point patterns (e.g., [26][27][28][29][30][31]). It is a statistical method developed by Silverman [47] and acts as a smoothing filter for discrete points (Equation (1)). The KDE method is based on the value of i at each sampling point (x, y) in a kernel. The density contribution of each sampling point was calculated using the kernel functions. The density value changes depending on the search radius range (bandwidth) and grid cell distance in each sample point. The method is already implemented in ArcGIS 10.2 from ESRI, Redlands, CA, U.S.A. and uses a quartic kernel, as proposed by Silverman [47]. In the search radius range, the sampling points closer to the grid cell center have a higher density contribution value. A cumulative density value for all kernels in the entire working area are computed using the method.
where; i = 1, . . . ,n are the input points. Only include points in the sum if they are within the radius distance of the (x, y) location. pop i : population field value of point I dist i : the distance between point i and the (x, y) location.
The parameter selection for KDE (e.g., kernel type and radius) was investigated in several studies (e.g., [48][49][50][51][52]) and yet there is no standard approach for this purpose. In this study, a spatial statistical analysis, average nearest neighbor (ANN) [53] implemented in ArcGIS 10.2, was used to analyze the distributions of the points and an approach based on these results was developed for the radius selection. The method computed the average nearest neighbor ratio and expected mean distances (D E ) and observed mean distances (D O ), as shown in Equation (2). where; DO is the observed mean distance between each feature and its nearest neighbor. DE is the expected mean distance for the features given in a random pattern. d i equals the distance between feature i and its nearest neighbor. n corresponds to the total number of features.
A is the area of a minimum enclosing rectangle around all features, or it is a user-specified area value [52].
If the ANN index is less than 1, the distribution is clustered. If the index is greater than 1, the points are dispersed. The results of the ANN analysis for each group are provided in Table 4. The flora points are clustered and the D E represents the mean sizes of the clusters. Due to the random distribution of bird and reptile points, the D O values are considered to be representative. On the other hand, there are only five observations (points) for mammals and they are dispersed. Since there is only one species type (Capra aegagrus Erxleben) in this group, its home range was obtained from [54] as 150-400 ha. It should be noted that identifying home ranges precisely is an exhaustive process and a simple and easy-to-use approach is needed for planning purposes. In this study, a 1000 m kernel radius was selected for all groups based on the D O and D E values, and the home range obtained from [54]. In addition, during the initial investigations on kernel radius, it was observed that using a smaller radius results in a disconnected and highly fragmented ecological network. On the other hand, using a much larger kernel radius would cause over-smoothing of the data. A 1000 m radius was found to be visually optimal and it also results in a ca. 3.14 ha home range for the mammals. The following strategy can be derived from the current implementation based on the ANN analysis: If distribution pattern is clustered → check D E If distribution pattern is random → check D O If distribution pattern is dispersed → check home range for fauna (if possible) or EUNIS habitat class for flora. The natural breaks classification is used to reclassify the output KDE values of each group into five categories. The method implemented in ArcGIS 10.2 uses the Jenks optimization algorithm [55] that minimizes the squared deviations of the class means. After this step, the classes were considered as habitat zones, as shown in Table 5. The four habitat classes were combined using an integrated density index approach to obtain a joint habitat pattern. The approach used here compares the four class values (flora, mammals, reptiles and birds) for every grid point and takes the maximum class value as output. Although other integration approaches (e.g., averaging the class values, calculating the sums of all class values for every pixel and reclassifying into five classes again, etc.) could be used; the proposed approach allows us to represent the importance level of species when having a small number of samples (e.g., mammals) better and therefore is preferred. The results of different bandwidth tests with 500 m and 800 m for endemic flora and a single class for fauna data (instead of three different classes for mammals, reptiles and birds) are presented in Figure 10. These radii cause disconnected and fragmented ecological network. On the other hand, as can be seen from Figure 10b, the core zones of threatened mammals (denoted with red triangles) cannot be represented in the final map when all fauna data are grouped together. The KDE analysis results of all four groups classified with natural breaks are provided in Figure 11. The same anomaly occurs if other integration approaches (e.g., averaging the class values, etc.) were applied and therefore using maximum class value allows the optimal representation of the core zones, buffer zones, etc., as presented in the following section. The natural breaks classification is used to reclassify the output KDE values of each group into five categories. The method implemented in ArcGIS 10.2 uses the Jenks optimization algorithm [55] that minimizes the squared deviations of the class means. After this step, the classes were considered as habitat zones, as shown in Table 5. The four habitat classes were combined using an integrated density index approach to obtain a joint habitat pattern. The approach used here compares the four class values (flora, mammals, reptiles and birds) for every grid point and takes the maximum class value as output. Although other integration approaches (e.g., averaging the class values, calculating the sums of all class values for every pixel and reclassifying into five classes again, etc.) could be used; the proposed approach allows us to represent the importance level of species when having a small number of samples (e.g., mammals) better and therefore is preferred. Table 5. Natural breaks classification of the kernel density estimation (KDE) patterns.

Class Value Kernel Density Estimation (KDE) Spatial Configuration
5 High density pattern Core habitats 4 Medium high-density pattern Buffer patches 3 Medium-density pattern Ecological corridors 2 Medium low-density pattern Stepping stones 1 Low density pattern Isolated area The results of different bandwidth tests with 500 m and 800 m for endemic flora and a single class for fauna data (instead of three different classes for mammals, reptiles and birds) are presented in Figure 10. These radii cause disconnected and fragmented ecological network. On the other hand, as can be seen from Figure 10b, the core zones of threatened mammals (denoted with red triangles) cannot be represented in the final map when all fauna data are grouped together. The KDE analysis results of all four groups classified with natural breaks are provided in Figure 11. The same anomaly occurs if other integration approaches (e.g., averaging the class values, etc.) were applied and therefore using maximum class value allows the optimal representation of the core zones, buffer

The PA Network
The integrated results clipped using the coastline and the PA boundary are presented in Figure  12. Unlike using discrete points for demonstrating the habitats, the results indicate a continuous pattern as expected, and show that the PA must be protected against habitat fragmentation and losses caused by it [56,57]. The five levels were interpreted as core habitat zones, buffer patches, ecological corridors, stepping zones and isolated elements for endemic and threatened species (Table 5). It should be noted that the collected data are still likely incomplete, due to the time span of the field campaigns, although the species were precisely identified and their locations were accurately measured [36]. Undersampling is an expected limitation of biodiversity determination projects due to the data collection methods (e.g., often manually by human interpretation) and the discrete nature of the measurements (discrete time intervals). The results mainly point out the areas for strict protection and the remaining parts should be analyzed by experts, if possible, and interpreted carefully by planners.

The PA Network
The integrated results clipped using the coastline and the PA boundary are presented in Figure  12. Unlike using discrete points for demonstrating the habitats, the results indicate a continuous pattern as expected, and show that the PA must be protected against habitat fragmentation and losses caused by it [56,57]. The five levels were interpreted as core habitat zones, buffer patches, ecological corridors, stepping zones and isolated elements for endemic and threatened species (Table 5). It should be noted that the collected data are still likely incomplete, due to the time span of the field campaigns, although the species were precisely identified and their locations were accurately measured [36]. Undersampling is an expected limitation of biodiversity determination projects due to the data collection methods (e.g., often manually by human interpretation) and the discrete nature of the measurements (discrete time intervals). The results mainly point out the areas for strict protection and the remaining parts should be analyzed by experts, if possible, and interpreted carefully by planners.

The PA Network
The integrated results clipped using the coastline and the PA boundary are presented in Figure 12. Unlike using discrete points for demonstrating the habitats, the results indicate a continuous pattern as expected, and show that the PA must be protected against habitat fragmentation and losses caused by it [56,57]. The five levels were interpreted as core habitat zones, buffer patches, ecological corridors, stepping zones and isolated elements for endemic and threatened species (Table 5). It should be noted that the collected data are still likely incomplete, due to the time span of the field campaigns, although the species were precisely identified and their locations were accurately measured [36]. Undersampling is an expected limitation of biodiversity determination projects due to the data collection methods (e.g., often manually by human interpretation) and the discrete nature of the measurements (discrete time intervals). The results mainly point out the areas for strict protection and the remaining parts should be analyzed by experts, if possible, and interpreted carefully by planners.

Interactions between the PA Network and the Human Use
The property ownership rights should be taken into account in a PA for effective protection. The private property areas are of particular interest, since protection should be carried out with the mutual agreement of the landowners. The protection-use balance can be achieved with the voluntary cooperation of the stakeholders [58][59][60][61][62] and the conservation strategies should be developed by social and financial analysis of the land property ownership and land use patterns. The distributions of different ownership types (i.e., private properties, private properties obtained after deforestation, forests under state ownership, and non-forest state properties), as shown in Figure 6, demonstrate the serious fragmentation of land registered under private ownership, both private parcels and newly registered parcels after deforestation. A major conservation strategy can be expropriation by the government in strict protection areas. The private ownership percentage (~20%) given in Table 6 also indicates the need for expropriation for nature conservation. For this purpose, a preliminary analysis was made using the privately owned parcels, the habitat network obtained from the developed methodology, and the EUNIS classes. The analysis extracts the privately owned parcels from the dataset that fall inside the core habitat zones and the buffer patches (red polygons in Figure 13). However, parcel areas that fall into the EUNIS classes of residential areas (villages), cemeteries, agricultural constructions and mixed crops of market gardens and horticulture were excluded from the analysis. According to the national regulations [33] for the use permissions in unplanned areas, the borders of these four classes were expanded by 100 meters prior to the analysis. Although the results given in Figure 13 show the areas that need to be expropriated, the remaining parcels should still be examined by experts for a final decision, due to the complex nature of the biodiversity data.

Interactions between the PA Network and the Human Use
The property ownership rights should be taken into account in a PA for effective protection. The private property areas are of particular interest, since protection should be carried out with the mutual agreement of the landowners. The protection-use balance can be achieved with the voluntary cooperation of the stakeholders [58][59][60][61][62] and the conservation strategies should be developed by social and financial analysis of the land property ownership and land use patterns. The distributions of different ownership types (i.e., private properties, private properties obtained after deforestation, forests under state ownership, and non-forest state properties), as shown in Figure 6, demonstrate the serious fragmentation of land registered under private ownership, both private parcels and newly registered parcels after deforestation. A major conservation strategy can be expropriation by the government in strict protection areas. The private ownership percentage (~20%) given in Table 6 also indicates the need for expropriation for nature conservation. For this purpose, a preliminary analysis was made using the privately owned parcels, the habitat network obtained from the developed methodology, and the EUNIS classes. The analysis extracts the privately owned parcels from the dataset that fall inside the core habitat zones and the buffer patches (red polygons in Figure 13). However, parcel areas that fall into the EUNIS classes of residential areas (villages), cemeteries, agricultural constructions and mixed crops of market gardens and horticulture were excluded from the analysis. According to the national regulations [33] for the use permissions in unplanned areas, the borders of these four classes were expanded by 100 m prior to the analysis. Although the results given in Figure 13 show the areas that need to be expropriated, the remaining parcels should still be examined by experts for a final decision, due to the complex nature of the biodiversity data.  Figure 13. Cadastral analysis results and the determined habitat network. Infrastructure facilities in rural areas provide economic development but may destroy natural assets [63] and threaten the biodiversity of PAs [64][65][66][67]. Roads affect the PAs by their width, traffic level and speed [68]. Amphibians, reptiles, birds and mammals are mostly affected by the roads [69,70]. The impact of roads on the ecological network is in the range of 100-600 m [65] and according to different road types (provincial, district and rural roads), the impact area can be determined as 500, 100 and 50 m, respectively [71][72]. The roads in the PA were buffered using different radii (500, 100 and 50 m) based on their type (province, district, rural) and visualized together with the determined ecological network for defining the conservation strategies and to support the decision making for the construction of new roads. Figure 14 shows the road impact zones (buffers) and settlement areas overlaid with the determined ecological network. From the Figure 14, it can be seen that they have impacts on the core habitats and the buffer patches. Infrastructure facilities in rural areas provide economic development but may destroy natural assets [63] and threaten the biodiversity of PAs [64][65][66][67]. Roads affect the PAs by their width, traffic level and speed [68]. Amphibians, reptiles, birds and mammals are mostly affected by the roads [69,70]. The impact of roads on the ecological network is in the range of 100-600 m [65] and according to different road types (provincial, district and rural roads), the impact area can be determined as 500, 100 and 50 m, respectively [71,72]. The roads in the PA were buffered using different radii (500, 100 and 50 m) based on their type (province, district, rural) and visualized together with the determined ecological network for defining the conservation strategies and to support the decision making for the construction of new roads. Figure 14 shows the road impact zones (buffers) and settlement areas overlaid with the determined ecological network. From the Figure 14, it can be seen that they have impacts on the core habitats and the buffer patches.

Conservation Strategies
The strategies reported by biodiversity researchers [36] mainly include observations on biodiversity elements and actual land use within different parts of the PA with qualitative

Conservation Strategies
The strategies reported by biodiversity researchers [36] mainly include observations on biodiversity elements and actual land use within different parts of the PA with qualitative geographical definitions (e.g., west of Kekova Island, west of the PA, etc.). Although these recommendations describe the biodiversity and land use status, the PA management strategies and regulations are mainly defined by the MoEU. The strategies and regulations can be summarized as: • A PA should be determined based on the assessments of ecological, geological, geomorphological, hydrological and landscape characteristics in the area.

•
The PA strategies, which should support the decision makers, must be based on qualitative and quantitative assessment results.

•
The biological/ecological importance of a PA should be exhibited with the existence of critical species and the actual status and pattern of their habitats. The main aim of biodiversity studies carried out under the direction of the MoEU is to determine the biodiversity (flora and fauna) and habitat characteristics. The ecological network should also be demonstrated with such studies.

•
The whole ecosystem should be assessed together with the spatiotemporal distributions of critical species and their habitats, and also their relations.

•
Thus, planning and protection can be based on ecological studies and scientific data.
Here, these principles were expanded for Kas-Kekova PA based on the outcomes of this study. After the determination of the ecological characteristics and the PA habitat network, conservation-use principles for the protection of endangered species and sustainable land use could be proposed. First, joint efforts were carried out with the government authorities responsible for the development of the Kas-Kekova Conservation Plans. For sustainable land use, the following conservation strategies can be proposed here:

•
The ecological network should be integrated into the land use plans. Strict protection measures, such as no construction, must be taken in the core habitat zones. Restoration of the core habitats should be improved.

•
Buffer patches (zones) are helpful to maintain the core zones and they contain important habitats. Buffer patches ensure distinctive habitat isolation and must be defined as strictly PA as well.

•
Infrastructure facilities in rural areas provide economic development, but easily lead to the destruction of natural assets [58][59][60][61][62]. To provide effective land management in the PA, the expansion of the settlements should be prevented.

•
Land consolidation should be carried out in the PA to support the protection of the habitat and obtain an optimal land use pattern. Existing orchards can be supported with ecological agriculture and used as connecting elements for the other zones, especially for endangered species and threatened habitats.

•
The Government should expropriate the private property areas that are located in the core habitats and the buffer patches. If the expropriation process cannot be executed rapidly, their owners can be supported by the social system to be able to prohibit any economic activity (e.g., tourism, farming, etc.) or the farmers can be trained and encouraged for ecological agriculture activities to protect and enhance the environment on their farmland.

•
The roads in the PA should be analyzed based on their type (province, district, rural) and compared with the determined ecological network for defining the conservation strategies and supporting the decision making for the construction of new roads and the reconstruction of existing ones.

Conclusions and Future Works
A quick and feasible method based on KDE and the integrated density index to define the habitat patterns, buffer zones and ecological corridors of threatened species in Kas-Kekova PA, Turkey was investigated in the present study. According to the results of spatial analysis and interpretation, a protection strategy was developed and proposed as a comprehensive approach for effective land management. From the results, it can be concluded that the KDE is powerful to derive a connected habitat network from discrete points of threatened species locations, although the selection of a kernel radius that may be suitable for diverse species types can be challenging. Here, ANN was employed to analyze the point patterns (random, clustered, disperse) and a kernel radius selection approach is proposed based on the results. However, the ANN method may not be suitable with few samples (e.g., for endangered mammals) and more expertise and field studies are needed for the acquisition of more knowledge about the threatened species and thus about the kernel radius selection. The integrated density index approach was used for combining the KDE results of different point groups (i.e., mammals, reptiles, birds and flora) and using the maximum input value for each output index value avoids over-smoothing and enables the representation of the importance value (e.g., core zones) of species with small numbers of samples. The output of the methodology is a base map which demonstrates the ecological network of threatened flora and fauna that can connect the fragmented habitats and be used to analyze the habitat patterns.
Nature conservation activities are often limited by the property ownership rights and the actual human use (e.g., roads) in PAs. As stated in the conservation strategies, expropriation and land consolidation tools should be developed by the local authority to improve this situation. The biodiversity studies must be integrated into the regional land use plans in Pas, since the land use plans are the most important means for sustainable nature conservation. For effective land management in a PA, a spatial environmental plan must be prepared by giving the priority to the biodiversity reports. If the priorities are not clarified with biodiversity maps, urban planners could easily get away from the protection purpose or cannot determine the land use boundaries in spatial environmental plans.
Despite the importance of biodiversity determination projects in PAs, they require a high level of expertise, extensive manual labor, and field works which should be spread over many seasons and years. On the other hand, with the recent advances in spatial data collection, interpretation and analysis methods-such as volunteer geographical information (VGI), citizen science (CitSci), remote sensing, new sensor technologies, and GIS as data storage, analysis, sharing and visualization platform-the task becomes easier. The increase in the studies on biodiversity and land use monitoring with remote sensing techniques indicates that PA planning and management will also be supported better with the help of geospatial technologies. Future work in Kas-Kekova PA and other PAs in Turkey should involve the development of a system-based monitoring approach complementary to the biodiversity determination projects.
Kas-Kekova PA is also a marine protected area. In this study, only the terrestrial habitats of the PA have been assessed. In future works, the analyses should be expanded to the marine part of the PA and a combined set of conservation strategies should be developed for the whole region, and the terrestrial and marine PAs should be evaluated together.