Natural Protected Areas as Providers of Ecological Connectivity in the Landscape: The Case of the Iberian Lynx

: The design of conservation plans for the improvement of habitats of threatened species constitutes one of the most plausible possibilities of intervention in the structure and composition of the landscape of a large territory. In this work we focus on the Iberian lynx in order to estab-lish potential ecological corridors using ecoinformatic tools from the GIS environment to improve connectivity between the existing natural spaces within the scope of its historical distribution. We processed 669 records of the presence of the lynx and six predictor variables linked to the habitat of the species. With this, corridors have been generated between natural areas. The determination of possible bottlenecks or dangerous areas (e.g., hitches on highways) allows for focusing efforts on their conservation. This type of approach seeks to improve efﬁciency in the design of measures aimed at expanding the territory’s capacity to host its populations, improving both its viability and that of all the other species that are linked to it. The proposals for action on the speciﬁc areas deﬁned by the models elaborated in this work would imply interventions on the land uses and existing vegetation types in order to improve connectivity throughout the territory and increase the resilience of its ecosystems.


Introduction
Human activities shape and transform territories on a global scale through the modification of the forms and properties of the surface and subsoil; being, therefore, a dominant factor in the evolution of the landscape in this current period [1,2]. As a consequence of this transformation, severe alterations in ecosystem dynamics are generated, resulting in population and species reductions and extinctions. In this context, the fragmentation of habitats and populations appears as one of the main factors responsible for the serious and rapid loss of biodiversity [3,4].
Habitat fragmentation results in a loss of connectivity in the landscape. The preservation and restoration of connectivity has become one of the main conservation objectives [5]. Understanding the ecological processes which depend on connectivity and taking effective planning measures requires an understanding of how landscape features can affect it. Taylor et al. [6] define connectivity as the level of facilitation or resistance to the movement of organisms between patches of habitat with resources. This resistance is determined by a landscape matrix which can present different levels of alterations and is able to modulate connectivity between patches or provide resources to species. In this way, a well-preserved matrix with a low degree of alteration can act as a buffer zone for habitat patches and mitigate the isolation of these. In addition, the matrix between these patches is characterized The main objective of this study is to provide a methodological approach to landscape assessment in order to maintain ecological connectivity, facilitating adequate landscape and land use management. We present this approach through the characteristics and habitat needs of a species of great interest in the Iberian Peninsula: the Iberian lynx. In this way, we aim to connect natural protected areas with good habitat and protection characteristics through a landscape matrix of very different qualities. To find out how the landscape influences the conservation of this species, we have set ourselves the following objectives: (a) to find out which factors most determine the distribution of the Iberian lynx in its historical range, (b) to generate a map of the potential distribution of the species in order to determine the suitability of the geographical area for its presence, (c) to evaluate different protected natural areas integrated into the historical distribution area of the species according to different factors: suitability of habitat and human influence, and (d) to evaluate the capacity for connection between these spaces and the generation of potential ecological corridors integrated into the landscape matrix that could be used by the species for its dispersion and conquest of future new territories.

Study Area
The study area covers approximately 176,254 km 2 located mainly in the central-west and south-west of the Iberian Peninsula ( Figure 1) where there is historical evidence of the presence of Iberian lynx. The landscape is very heterogeneous and combines mountainous zones belonging to the Central System, Montes de Toledo, Sierra Morena, and a large part of the Sub-Baetic Systems with different river basins (Guadalquivir, Guadiana and Tajo), especially the depression formed by the river Guadalquivir.
Sustainability 2021, 13, x FOR PEER REVIEW 3 of 17 that are randomly situated within the landscape are located in accessible habitat areas in a set formed by habitat patches and the links (connections) between them [36].
The main objective of this study is to provide a methodological approach to landscape assessment in order to maintain ecological connectivity, facilitating adequate landscape and land use management. We present this approach through the characteristics and habitat needs of a species of great interest in the Iberian Peninsula: the Iberian lynx. In this way, we aim to connect natural protected areas with good habitat and protection characteristics through a landscape matrix of very different qualities. To find out how the landscape influences the conservation of this species, we have set ourselves the following objectives: (a) to find out which factors most determine the distribution of the Iberian lynx in its historical range, (b) to generate a map of the potential distribution of the species in order to determine the suitability of the geographical area for its presence, (c) to evaluate different protected natural areas integrated into the historical distribution area of the species according to different factors: suitability of habitat and human influence, and (d) to evaluate the capacity for connection between these spaces and the generation of potential ecological corridors integrated into the landscape matrix that could be used by the species for its dispersion and conquest of future new territories.

Study Area
The study area covers approximately 176,254 km 2 located mainly in the central-west and south-west of the Iberian Peninsula ( Figure 1) where there is historical evidence of the presence of Iberian lynx. The landscape is very heterogeneous and combines mountainous zones belonging to the Central System, Montes de Toledo, Sierra Morena, and a large part of the Sub-Baetic Systems with different river basins (Guadalquivir, Guadiana and Tajo), especially the depression formed by the river Guadalquivir.

Study Species
The Iberian lynx is a slender cat, approximately 1 m long and 8-15 kg in weight [37]. It is an emblematic Iberian species considered to be one of the most seriously threatened cat species in the world and is currently listed as "endangered" according to the International Union for Conservation of Nature (IUCN) categories [38,39]. The Iberian lynx is

Study Species
The Iberian lynx is a slender cat, approximately 1 m long and 8-15 kg in weight [37]. It is an emblematic Iberian species considered to be one of the most seriously threatened cat species in the world and is currently listed as "endangered" according to the International Union for Conservation of Nature (IUCN) categories [38,39]. The Iberian lynx is considered a strict species in terms of habitat requirement, being found in areas between 400 and 1300 m above sea level and closely linked to the Mediterranean mountain range, especially the extensive and dense scrub [40]. It has an elusive nature, avoiding open areas with human influence in all stages of its biological cycle [41]. In their dispersal phase, lynx are able to use remnants of suitable habitat, formed mainly by Mediterranean scrub, as "steps" to travel through a fragmented landscape [42] and through a poorer quality matrix formed by woodlands or forest plantations of pines or eucalyptus [43].
The reduction in its range has been linked to several key factors: on the one hand, agricultural and forestry transformations; on the other hand, the decline of its main prey, the rabbit (Oryctolagus cuniculus Linnaeus, 1758) due to myxomatosis and viral hemorrhagic disease [40,41,44]. Finally, the increase in non-natural mortality of this species caused mainly by illegal hunting and road traffic accidents is a serious obstacle to its recovery [40,42,45].

Species Data
The Iberian lynx presence (n = 900) was obtained from the Global Biodiversity Information Facility (GBIF: https://www.gbif.org) [46] through human observation. Most of these records are integrated in the dataset of the National Biodiversity Inventory (2007) of the Ministry of Environment, Rural and Marine Affairs, collected between 1980 and 2007. The data were processed to eliminate duplicates and incorrect records, collecting a total of 669 presence records at 10 × 10 km resolution.

Environmental Data
Twenty-two variables were first collected relating to climatic, topographical, and land use characteristics. The trophic variable (rabbit density) was not used due to the lack of a thematic cartography in this study.
The bioclimatic variables were obtained from WorldClim-Global Climate Data version 2.0: http://worldclim.org/version2 for the period 1970-2000 at a resolution of 1 × 1 km. The slope was derived from a Digital Elevation Model (DEM) at 200 × 200 m resolution from the IGN-National Geographic Institute: https://www.ign.es, using the ArcGis Slope function [47]. The vegetation cover variable was obtained from the Spain Forest Map, project carried out during the years 1997 to 2006 and available at Ministerio para la Transición Ecológica (MITECO Spanish acronym, https://www.miteco.gob.es). We extracted the vegetation cover using raster calculator (ESRI, 2016).
Six of those 22 variables were finally selected following bibliographic sources on the biology and ecology of the Iberian lynx [41,45]. They are related to climatic, topographical, and land use aspects (Table 1). Seasonal precipitation (variation coefficient) [50] All environmental variables were standardized at geographic coordinates (Datum WGS-1984) and resampled at a spatial resolution of 10 × 10 km using bilinear interpolation sampling [47]. Previously, all continuous variables were checked for collinearity effects, through the "ggpairs" function of the "GGally" package [51] in R version 3.1.3 [52]. Those variables with (|r| ≥ 0.75) were removed from the final set of predictor variables (Table S1).

Modelling Process with MaxEnt
Habitat suitability of Iberian lynx was modelled using MaxEnt version 3.4.1, a machinelearning process that uses presence-only data [25]. MaxEnt gives insight about what predictors are important and estimates the relative suitability of one place vs. another, as well as the probability of occurrence [53,54]. This approach has been demonstrated to perform well in a diverse set of modelling scenarios and is widely used in a great number of studies in ecology, biogeography and conservation [55][56][57][58].
Ten replicates were performed by Bootstrap procedure in "cloglog" format using 80% of the presence records as model training and the remaining 20% to evaluate the models. The background was of 10,000 random points from study area. Replicas used herein were performed using the MaxEnt default parameters [25]. Variable importance was calculated to assess the relevance of each predictor in the models and response curves were calculated to interrogate the relationship between the response (i.e., Iberian lynx presence) and each explanatory variable The models were evaluated using the area under the curve of the receiver operating characteristic (AUC). AUC measures the ability of a model to discriminate between sites with occurrence, versus those where it is absent (herein, background points). The AUC ranges from 0 to 1 (0.5 = random, 1 = perfect). Finally, the model was projected into the study area to generate a habitat suitability map of Iberian lynx. The criterion of maximizing the sum of sensitivity and specificity (maximum test sensitivity plus specificity) was used as a threshold to define the suitability of the cells [59,60].

Selection of Core Areas
The Natural Protected Areas (NPAs) in the study area have been used as core areas. In particular, national parks and natural parks with protection category that have associated regulations and conservation measures for species and habitats were included.
To improve the analysis, we decided to merge those geographically adjacent NPAs and consider them as a unique unit ( Figure S1). This is the case of Sierra Norte Sevilla, Sierra Hornachuelos, and Sierra de Aracena y Picos de Aroche (unit 2); Valle de Alcudia and Sierra Madrona, Sierra de Andújar, and Sierra de Cardeña y Montoro, (unit 4) and Quilamas and the Batuecas-Sierra Francia (unit 17). The selected NPAs (25) are shown in Table 2, where the numerical coding adopted is indicated in order to facilitate the interpretation and discussion of the results.

Resistance Raster
The next step in the connectivity analysis is the generation of a raster showing the resistance or suitability of the terrain-landscape matrix-for the movement of the species [30,61]. A specific ArcGis extension called Corridor Designer (http://corridordesign.org) was used for this purpose [62], which allows for the design of ecological corridors based on a resistance raster.
This raster requires a previous weighting by the researcher of different selected variables, in this case: altimetric information representing the relief in the area (DEM) [48], topographic slope (Slope), a set of global data covering human impact (population density, land use, infrastructure, etc.) represented by the Human Influence Index (HII) was downloaded from the SEDAC-Socioeconomic Data and Applications Center database: https://sedac.ciesin.columbia.edu [63] and land cover information obtained from the National Forest Inventory through the Forest Map [49].

Corridor Design
The resistance map is used by the Corridor Designer as a basis for calculation of the corridors with a range of 0 to 100, where 0 represents the most restrictive value for the species and 100 the optimum. The resolution of this habitat model was 1 × 1 km. The lowest cost corridors are generated through a cost-distance analysis. Therefore, they are distance based corridors that minimize the cumulative cost of each start cell location through the above mentioned resistance raster. In this way the cost of making a move through each cell is defined. Table 2. Natural Protected Areas (NPAs) selected as core areas to be connected in the corridor network. An alphanumeric coding is shown for a better understanding of the maps, the Protection Figure (Figure NPA), its area in km 2 , the average suitability, the Human Influence Index (HII) and information about whether they currently host (Yes or No) an Iberian lynx population. CP: Current Population. NP: National Park; NatP: Natural Park; RegP: Regional Park; ZRI: Zone of Regional Interest; SCZ: Special Conservation Zone. Additionally, Corridor Designer considers biological information of the species through the following parameters: the threshold value, the minimum breeding area and the minimum host area of a population. The threshold value is a suitability value that discriminates between habitats where reproduction is possible and those where it is not [62]. It is related to the characteristics of the species and its habitat requirements. Majka et al. [62] recommend a standard threshold value of 60, but in our case, we selected a more restrictive value of 80 due to the narrow ecological requirements of the species. The minimum breeding area corresponds to a surface value large enough for the species to reproduce. The minimum breeding area for a population is defined as an area large enough to support reproduction for 10 years or more. The minimum breeding area and the minimum hosting area of a population adopt values of 500 and 2500 hectares, respectively (based on the connectivity study carried out by Puerto Marchena and Muñoz Reinoso [20]). At this stage of the process, this information is required to calculate suitable patches within the areas to be connected and use them as starting and finishing points, prioritizing the population patches.

Cod
Depending on the surface of the landscape considered, nested polygons (corridors) with different surfaces (polygon-corridor 1%, 2%, 3%, . . . ) are generated as possible alternatives for transit between the nodes to be connected and can be considered as buffer areas. These are low-cost corridors that widen in quality areas and narrow where the habitat is unsuitable [64].

Corridor Evaluation
Once the ecological corridors are defined, we attempted to evaluate and identify those landscape elements that could be most vulnerable in order to prioritize their conservation or restoration. The identification of these elements or critical areas is summarized as follows: (1) analysis of bottlenecks along the corridors using the Bottleneck Analysis tool, (2) characterization of the core areas according to the index of human influence and the predicted average suitability, and (3) determination of the contribution of the landscape elements to the overall connectivity of the network using the Conefor sensinode 2.6 software.
Bottleneck analysis was performed with the Bottleneck Analysis tool of the Corridor Designer Evaluation extension [62]. This tool allows the location of bottlenecks along corridors, which are narrows that can make the corridor less effective. This is important for those species sensitive to edge effect, infrastructure, etc. For this purpose, we tested with a threshold value of 1500 m.
In addition, the core areas and corridors were classified according to habitat suitability (Maxent's model) and Human Influence Index (HII). On the other hand, we considered the habitat patches generated with Corridor Designer. The existence of quality habitat patches along the corridor may facilitate the transit of the species. This is fundamental in those areas where the corridor crosses areas of low-quality habitat [42], and also for those long corridors functioning as stepping stones. This could serve to facilitate measures and proposals for new reserves or protected areas in a strategic way.

Prioritization of Landscape Elements for Conservation and Restoration
The importance of landscape elements to the connectivity of protected areas for the Iberian lynx is quantified using a connectivity metrics approach, in our case, the Probability of Connectivity (PC) [36]. We used Conefor sensinode 2.6 software [65] available and updated at www.conefor.org.
Two input files are required: one relating to the nodes and the other to the links. In the first (nodes), the different patches to be connected are collected and each one of them is assigned an attribute, which is usually the area. We used the natural break method or Jenks optimization method [66] to classify in three intervals, both the size and the average suitability of the core areas. This method reduces the variation of each class, while maximizing the variance between classes. These areas were valued from 1 to 3 depending on the interval in which they have been located. Finally, we summed these values, obtaining values between 2 and 6. Thus, we obtained an attribute that integrates the surface and the average suitability (obtained through the species distribution model). The link file includes the references of the nodes that are connected and the length between them. This length is usually Euclidean, but we used the minimum cost distances obtained from the resistance map.
The PC approach is often used in connectivity analyses [67][68][69]. We focus on the negative effects for dispersion and re-colonization resulting from the loss or degradation of landscape elements. The PC index can be defined as the probability that two points (animals) that are randomly situated within the landscape are located in accessible habitat areas in a set formed by habitat patches and the links (connections) between them. The landscape elements are classified according to their contribution to the overall availability and connectivity of the habitat. This contribution is calculated by the percentage of variation of the PC index obtained after the individual elimination of each element [36].
For this index it is possible to differentiate three separate fractions that quantify the ways in which elements contribute to the overall connectivity and habitat availability of the landscape: dPCintra(k) measures intrapatch connectivity, while dPCconnector(k) and dPCflux(k) measure connectivity between patches in relation to a given landscape element k. Links only contribute to connectivity through the dPCconnector(k) fraction. Therefore, a patch will be more or less important (dPCk) depending on the intrinsic characteristics and its topological position within the network [33,70] Finally, the calculation method of dPC k and dPC k connectors was based on the value of the mean, the median or maximum dispersion distance. We tested with a value of 60 km so that the dispersion probability of an individual travelling 60 km in his dispersive phase will be 0.5. Another interesting metric is the Betweenness Centrality metric (BC), which refers to the degree to which the optimal patch paths pass through a particular node k. It provides a more complete view of how patches can be providers of connectivity [33,70].
We also carried out an analysis with Conefor sensinode to determine the elements of the proposed corridor network that contribute most to the availability and connectivity of the habitat. This was analyzed with the dPCconnector fraction. We selected the priority elements with the highest values of overall importance, and represented the results in ArcGis 10.5.

Evaluation of the Suitability Model
The final SDM showed AUC values of 0.779 ± 0.014 (average ± standard deviation for 10 replicates of the models). An AUC > 0.7 constitutes a result with good discriminatory power [71], indicating that the model is better than random, and that the suitable environmental areas predicted by MaxEnt are highly correlated with presence records. The method for selecting the suitability threshold is a very important step in the final processing of SDMs and depends fundamentally on the objective of the study. Liu et al. [60] consider that if the purpose is the search for new populations of a species, it will be of greater importance to minimize errors of commission (false positives), while if the objective is the establishment of a conservation area, it is possible to tolerate greater errors of commission. Thus, in our work the selection of suitable areas has been carried out according to the threshold established through the criterion: maximum test sensitivity plus specificity, with a value of 0.416 (41.6%) and two suitability maps have been drawn up (Figure 2), one in continuous and the other in binary (suitable-not suitable).  The variables that best explain the distribution of the lynx are, in this order: slope, vegetation cover, and the seasonal precipitation (Bio15), followed by the seasonal temperature (Bio4) and the annual precipitation (Bio12), and with a lower contribution, the average annual temperature (Bio1). The environmental variables with the greatest gain when used in isolation are slope, vegetation cover, and Bio15 ( Figure 3) because they provide the most useful information on their own. The variables that best explain the distribution of the lynx are, in this order: slope, vegetation cover, and the seasonal precipitation (Bio15), followed by the seasonal temperature (Bio4) and the annual precipitation (Bio12), and with a lower contribution, the average annual temperature (Bio1). The environmental variables with the greatest gain when used in isolation are slope, vegetation cover, and Bio15 ( Figure 3) because they provide the most useful information on their own.
As can be shown on the map of MaxEnt (Figure 2), two large areas with high suitability values can be seen: one that crosses the study area from SE to NW and another, more isolated, in the south of the study area. On the other hand, there is a large area of low suitability that mainly affects the west central part of the area under analysis and is related to the large presence of large crop cultivation areas on the banks of the Guadiana river and Tierra de Barros [72]. The case of Doñana is similar: connectivity with the western Sierra Morena is limited as it is located in a hostile environment, with extensive areas of  The variables that best explain the distribution of the lynx are, in this order: slope, vegetation cover, and the seasonal precipitation (Bio15), followed by the seasonal temperature (Bio4) and the annual precipitation (Bio12), and with a lower contribution, the average annual temperature (Bio1). The environmental variables with the greatest gain when used in isolation are slope, vegetation cover, and Bio15 ( Figure 3) because they provide the most useful information on their own. As can be shown on the map of MaxEnt (Figure 2), two large areas with high suitability values can be seen: one that crosses the study area from SE to NW and another, more isolated, in the south of the study area. On the other hand, there is a large area of low suitability that mainly affects the west central part of the area under analysis and is related to the large presence of large crop cultivation areas on the banks of the Guadiana river and Tierra de Barros [72]. The case of Doñana is similar: connectivity with the western Sierra Morena is limited as it is located in a hostile environment, with extensive areas of There are areas of high suitability along the Central System, from the Sierra de Gata to the Sierra de Gredos and the Tiétar Valley in the northwest of the study area. The lack of information in the Portuguese territory limits the analysis of the suitability in this geographical area, although there is a glimpse of an area that apparently would present high suitability along the Portuguese border, from the west central corner of the study area to the Sierra de Gata and the Rebollar in the province of Salamanca, passing through Serra Malcata (Portugal). Similarly, there is an area of high suitability along the Portuguese border in the southwest corner of the area under analysis, which we suppose could extend to areas such as the Vale do Guadiana (current lynx population) and the Algarve, both in Portugal. The human influence rate along this entire border strip is low, mainly in the northwest of the territory studied, so it can be said that in absence of information about rabbit density, it is an area that holds good conditions to host the lynx. In fact, it was considered by some authors [17,[73][74][75] to be an area of great importance and the last strongholds of the lynx population in the Iberian Peninsula.

Suitability of Natural Protected Areas
The 25 NPAs selected (Table 2) have an average surface area of 674 km 2 . The least extensive NPA is the Despeñaperros Natural Park with 76.5 km 2 , while the most extensive is Sierra de Cazorla, Segura, and Las Villas with 2100 km 2 .
The average suitability of the NPAs (between 0 and 100%) is 59% according to the suitability map obtained with the MaxEnt prediction. The most suitable NPAs are Sierra de Castril with 86%, Despeñaperros with 82%, and Sierra de Andújar with 82%. On the other hand, there is only one NPA that is less suitable than the threshold value, Embalse de Orellana and Sierra de Pela, with only 31%, whose inclusion in this study was due only to its strategic geographical location. The core areas have been represented on a map for connectivity analysis and were classified according to their suitability ( Figure 4 and Table 2).
The Human Influence Index (HII) [63] (WCS and CIESIN, 2005) has been evaluated independently for each NPA ( Table 2). The core areas present values of around 30%. The natural park of the Sierras Subbéticas presents the highest index of human activity, with values of 40.4%. On the other hand, Cabañeros National Park presents the lowest index, with 18.6%. The Human Influence Index (HII) [63] (WCS and CIESIN, 2005) has been evaluated independently for each NPA ( Table 2). The core areas present values of around 30%. The natural park of the Sierras Subbéticas presents the highest index of human activity, with values of 40.4%. On the other hand, Cabañeros National Park presents the lowest index, with 18.6%. Table 2. Natural Protected Areas (NPAs) selected as core areas to be connected in the corridor network. An alphanumeric coding is shown for a better understanding of the maps, the Protection Figure (Figure NPA), its area in km 2 , the average suitability, the Human Influence Index (HII) and information about whether they currently host (Yes or No) an Iberian

Corridors Evaluation
Corridor 2-3 has important bottlenecks, with 73.1% narrowing below the established threshold (1500 m). On average, 16.4% of the kilometers that run through the proposed corridors are in areas with widths of less than 1500 m (Table S2).
The corridors were classified according to the habitat suitability of the territory predicted by MaxEnt (Figure 4). The average suitability of all the corridors is 55.3%, with those with the highest suitability being 4-14, 17-16, and 6-5 with values above 70%. On the other hand, the three corridors with suitability levels below 40% are 13-12, 2-3, and 2-11.
Landscape elements were classified in three categories by using the natural breaks of Jenks. Only four nodes (NPAs)-Sierra Morena West (2), Sierra Morena East (4), Sierras Subbeticas (3) and Sierra de Cazorla (6) with values above 3.5%, highlighting Sierra Morena East with 9.7% and eight connectors (4-19, 6-7, 16-17, 1-2 and 3-9) with dPC values above 5%-were assigned to the superior category ( Figure 6). The average length (Table S2) of the 29 corridors is 79.2 km. Corridor 13-12 is the longest at 168.5 km and Corridor 17-18 is the shortest at 11.5 km. The average width of the 29 corridors is 3.89 km (calculated on the 1% polygon). The corridor with the smallest average width is Corridor 2-3 with an average width of 1.2 km. Corridors 3-10 and 6-19 have the greatest narrowing, with widths of 0.164 km. The average minimum width of the 29 corridors is around 1.7 km.

Corridors Evaluation
Corridor 2-3 has important bottlenecks, with 73.1% narrowing below the established threshold (1500 m). On average, 16.4% of the kilometers that run through the proposed corridors are in areas with widths of less than 1500 m (Table S2).
The corridors were classified according to the habitat suitability of the territory predicted by MaxEnt (Figure 4). The average suitability of all the corridors is 55.3%, with those with the highest suitability being 4-14, 17-16, and 6-5 with values above 70%. On the other hand, the three corridors with suitability levels below 40% are 13-12, 2-3, and 2-11.

Discussion
Our study is based on the use of existing protected natural areas included in the landscape matrix to generate a network of corridors that will allow the identification of natural areas, connectors, and zones where conservation and restoration work can be most interesting in terms of the conservation of the Iberian lynx and its potential dispersal. The latter is of great interest, since it will allow, if necessary, for greater ease on the part of administrations to carry out restoration and conservation measures. Our approach is complementary to others, particularly the interesting study carried out by Blazquez-Cabrera et al. [21], in which they analyzed the importance and priority of restoring ecological corridors established between current Iberian lynx populations in the Iberian Peninsula.

Discussion
Our study is based on the use of existing protected natural areas included in the landscape matrix to generate a network of corridors that will allow the identification of natural areas, connectors, and zones where conservation and restoration work can be most interesting in terms of the conservation of the Iberian lynx and its potential dispersal. The latter is of great interest, since it will allow, if necessary, for greater ease on the part of administrations to carry out restoration and conservation measures. Our approach is complementary to others, particularly the interesting study carried out by Blazquez-Cabrera et al. [21], in which they analyzed the importance and priority of restoring ecological corridors established between current Iberian lynx populations in the Iberian Peninsula.
With this type of study, we intend to combine the characteristics-fundamentally abiotic-of the landscape with the biological aspects of the species to promote their conservation.

Connectivity Studies in Environmental Assessments. The Practical Approach
Connectivity studies are essential in environmental impact assessments and strategic environmental assessments. They make it possible to discriminate between different alternatives in the design of large infrastructures that are promoted in the landscape, mainly communication routes (highways, roads, railways, etc.) [76]. In addition, they allow the characterization of possible wildlife passage areas to reduce barrier effects, as well as the identification of "black points" of wildlife mortality due to traffic accidents [77,78] and of natural habitat restoration areas that can be used as "steps" by wildlife. The graphic approach allows for the simplification of a landscape mosaic and its complex network of functional connections into elements that form nodes and links. The nodes represent habitat sites surrounded by hostile habitats and the links symbolize the dispersal capacity of a species between two nodes. This method allows us to represent the landscape as a network of interconnected patches and to perform complex computational analyses on the connectivity of the landscape [29,34,79]. Usually, landscape connectivity studies have a strong local character; however, nowadays these analyses are increasingly applied to large territories and regions [11]. These regional connectivity studies provide great value as support tools for the future management of the landscape matrix. There are interesting examples of large-scale connectivity studies such as the initiative known as Yellowstone to Yukon (Y2Y) launched in 1998 [80], or the analysis of the connectivity of forest ecosystems integrated into the Natura 2000 Network in Spain, promoted by de la Fuente et al. [68].
The results obtained from this proposal allow the identification of key landscape elements according to their importance, both for their intrinsic characteristics (human influence index, average suitability) and for their contribution to connectivity, either by analyzing bottlenecks or by quantifying their contribution to overall connectivity.

Quality of Natural Protected Areas
As already mentioned, the NPAs included in this study were selected for their adequate ecological integrity, their strategic geographical location, and their protection category, which allows for the implementation of conservation measures. This is important, as some studies affirm that protected areas have lower rates of habitat loss [81,82]. In our case, some of the NPAs that currently have a lynx population are the most protected, such as Doñana, Sierra de Andújar, or Valle de Alcudia. In this sense, it would be expected that the variable of human influence (HII) was linked to the protection level, and therefore to the distribution of lynx, as argued by Palomares et al. [41]. HII values are between 18.6% (Cabañeros National Park) and 40.4% (Sierras Subbéticas Natural Park). In general, the most suitable areas for lynx have low HII (less than 30%), while the most anthropized NPAs (high Hll) were the Sierras Subbéticas, Sierra Grande de Hornachos, Embalse de Orellana, and-surprisingly-Doñana, due to the occupation of human infrastructures and agricultural plots.

Connectivity between NPAs. The Important Role of Sierra Morena
Although further research is still needed on the dispersal behavior of the Iberian lynx, some studies suggest distances of approximately 30-40 km in its wandering [41,42], but there are individuals who may well exceed these values. The average length of our corridors is 79.2 km, although they are very different in size. Each of these corridors should not be interpreted as a continuous unit of territory with optimal conditions for the species, but as a route connecting patches of suitable habitat (functioning as stepping-stones). As suggested by Ferreras [42] lynx may use these remnants of suitable habitat as "steps" to travel through a fragmented landscape. Furthermore, in terms of the distance between the fragments, the fragmentation of a territory seems to be directly related to the dispersal capacity of its species. Thus, in Doñana, moderate fragmentation effects seem to incite dispersers to risk longer journeys, exploring larger areas [42]. This has also been previously verified in experimental studies with wood mice and capercaillies [83,84] where the rate of movement between patches increased with habitat fragmentation.
According to the graphical approach, in our study area, node four (Sierra Morena Oriental) is the one that contributes most to the global connectivity of the network. In this case, the loss of this NPA would imply the division of the network into two large parts that would result in a great loss of connectivity. Without this node, NPAs such as two, three and six, as well as the adjacent ones, would be totally disconnected from those located north of the Sierra Morena Oriental. The Sierras Subbéticas (3) are a clear example of how a core area can function as a stepping stone to connectivity, since its loss would generate important disconnections between the NPAs of the Baetic systems (9,8,10), and those of the western and eastern Sierra Morena (two and four). Connectors such as 1-2 and 2-3 contribute greatly to network connectivity. The first one does so mainly because its elimination would determine the disconnection of an important NPA such as Doñana National Park. Similarly, Corridor 2-3 is a vital connector, establishing the connection between the Western Sierra Morena and the Baetic and sub-Baetic systems. This would be in line with the results obtained by Blazquez-Cabrera et al. [21] in which they identify two main points of action where restoration efforts should be focused. This is the axis formed by the populations of Doñana, Sierra Norte, and Matachel and the axis formed by the populations of Andújar, Guarrizas, and Sierra Morena Oriental, which coincides to a large extent with the results we have obtained on priority of areas and connectors.
As mentioned above, the Sierra Morena oriental patch (4) is a strategic node that is vital for maintaining connectivity in the corridor network. This node acts as a link between two large components in the study area and without it, global connectivity would be severely reduced and divided. It has the highest values of PC index for the three fractions (intra, flow, and connector) and also of BC (PC) centrality, showing its great importance both for the maintenance of connectivity and for the availability of its habitat and its possible function as a source of dispersal flow. If we consider the other characteristics used for the categorization of NPAs used in this study, the three natural parks forming this node are characterized by low values of the human influence index (less than 30%) and high habitat suitability (around 70%). All of this shows the fundamental role of this node in maintaining lynx populations and its function as a source of individuals and genetic structure for other areas in the process of recolonization.

The Importance of Corridor Width
Moreover, narrowing along the corridors can limit the efficiency of the course of the species, generating bottlenecks. It is important to consider the number of bottlenecks that a corridor has in its path, as well as the proximity between them. If bottlenecks are not too long and are spaced far enough apart, they may not be considered barriers to movement.
In general, the corridors generated in our study area have bottlenecks in low proportions along the course. Only the corridor linking the Sierra de Hornachuelos (2) with the Sierras Subbéticas (3) stands out, as it presents important narrowing in a large proportion (more than 70%). This could be due to the poor quality of the landscape matrix for the Iberian lynx, with high levels of human influence due to the proximity of large population centers, such as the city of Cordoba, as well as different highways and land dedicated to dry farming.

Conclusions
The use of tools in the field of eco-informatics, such as GIS in combination with SDMs and connectivity analyses, is a strategy for planning and managing biodiversity at broad territorial scales. This study presents a methodological approach aimed at facilitating decision-making in the management of the biotic and anthropic components of a landscape, in order to favor the conservation of its biodiversity. Our contribution, focused on the Iberian lynx and oriented to the theoretical design of corridors that connect different natural areas, aims to contribute to providing biological coherence to future anthropic intervention plans in the southwest of the peninsula, a territory of great natural wealth that is subject to numerous threats derived from the current development model.
Our results show that there are large areas and natural protected areas that maintain good characteristics to be able to host the Iberian lynx, both in terms of environmental suitability and human influence. Many of these areas do not host current populations of the species, although they did in the past, such is the case in El Rebollar, Sierra de Francia or Monfragüe; places where reintroduction and habitat management strategies can be carried out to facilitate their presence. Similarly, the analysis of connectivity in the study area provides relevant information on the contribution of connectors and NPAs to the general connectivity of the landscape, detecting critical points for the flow of individuals between