Global patterns and potential drivers of human settlements within protected areas

Protected areas (PAs) represent one of our most important conservation strategies for halting biodiversity loss. The number of PAs has increased remarkably over the last few decades. Yet, biodiversity is still being lost at alarming rates, even within many of those PAs. Understanding the factors that influence the levels of human pressure within PAs remains a key objective. In this study, we examined the factors associated with the human settlements’ levels within the world’s PAs. Using the random forests technique, an ensemble machine learning method, and a vast number of PAs (81 100–137 523), we assessed the importance of nine factors, including the PAs’ management objective as reflected by their International Union for Conservation of Nature (IUCN) Category. The IUCN classifies PAs into six categories ranging from strict nature reserves to areas in which multiple human uses are permitted. The prevalent but untested assumption is that human settlements’ levels within PAs vary according to their management objective, with less strict PAs having higher levels. Our results, however, show that the differences between the categories were for the most part minor. The most important predictor of human settlements was accessibility measured as the time required to reach the PA from the nearest major city. These findings were consistent across all of the world’s subregions. Other less important factors included the extent of croplands within PAs, elevation, and slope. Our findings suggest that PAs nearer urban centers tend to have higher human settlements’ levels regardless of their other characteristics, such as management objective and year of establishment. Managing those PAs successfully will be necessary to achieve the post-2020 global biodiversity targets and will require conservation strategies that acknowledge and engage the local communities.


Introduction
Protected areas (PAs) are often described as the cornerstone of our global conservation efforts to halt biodiversity loss [1,2]. Their importance is reflected in their remarkable increase in numbers over the last few decades (figure 1). Presently, there are more than 200 000 PAs across the globe [3], covering approximately 15% of its terrestrial area and 7.8% of its oceans [4]. The number of PAs worldwide is expected to increase further over the next years. The signatory parties to the Convention of Biological Diversity are expected to finalize the post-2020 global biodiversity targets at the 15th Conference of the Parties (COP-15), scheduled to be held in China in October 2021. According to the updated zero-draft of the Post-2020 Global Biodiversity Framework, released last August [5], the new target for the PAs and other effective areabased conservation measures will likely be set at 30% of the planet, to be achieved by 2030.
Despite the global commitment in PAs and their remarkable worldwide increase (figure 1), biodiversity is still being lost at alarming rates [6], even within many of those PAs [7]. This is because not all PAs are effective in eliminating human pressure [8][9][10]. Recent global analyses have shown that one-third of the world's protected land is under intense human pressure [8] and that many of the PAs worldwide continue to experience increasing levels of human disturbance [9][10][11]. Hence, understanding the factors associated with human pressure within PAs remains a crucial topic in the environmental sciences and has important implications for achieving the post-2020 global biodiversity targets [5] and several of the United Nations (UN) Sustainable Development Goals [12].
Not surprisingly, there have been several studies exploring the potential drivers of human pressure within PAs [12,13]. A key factor mentioned often in the literature is the PA's management objective (box 1) as reflected by its IUCN Category [14]. The IUCN classifies PAs into six categories, ranging from strict nature reserves (i.e. category Ia) to areas in which multiple human uses are allowed (e.g. categories V and VI; box 1). The categories were designed to reflect a gradient of naturalness [15] and permissible human uses [16]. It is often assumed that less strict PAs, such as those in categories V and VI, are less effective in mitigating human pressure [17,18] because they allow by design higher levels of human presence and activities [19,20]. This assumption has resulted in a long debate [21] regarding the relative effectiveness and the conservation value of those PAs [17,[19][20][21].
Understanding whether the IUCN categories play a critical role in influencing the levels of human pressure within PAs remains an important research objective for the following crucial reason: governments worldwide have been showing a clear preference for PAs with fewer restrictions (figure 1). Such PAs are likely easier to establish and less controversial, since they impose fewer limitations on the local peoples [22] and since they do not necessarily involve the expropriation of privately-owned land. Moreover, less strict PAs are thought to address many of the socioeconomic injustices associated with stricter PAs, which have been shown, in many cases, to disrupt severely the livelihoods of the local communities [23,24]. As a result, the six categories have increased disproportionately over the last few decades, with strictly PAs (such as those in categories I-II) increasing the least and the rest of the PAs (including those with no IUCN category) increasing the most (figure 1). If the critics of the less strict categories are correct in their assertions [17,19,20]-and if the levels of human pressure indeed vary across the categories-then the substantial increases in PAs worldwide may not translate into significant biodiversity benefits without appropriate and socially equitable conservation strategies, which account for this pattern [25,26].
Currently, the results in the literature regarding the categories remain inconclusive [27]. While some studies have found that strict PAs are more effective at mitigating human pressure [28,29], others have found the opposite [30,31], and yet others that there are no clear differences between the categories [2,32]. The majority of these studies, however, have focused on different geographic regions and have used dissimilar analytical methods and measures of human pressure [27]. Consequently, their findings are not always comparable [27]. A recent global study [27], based on more than 19 000 PAs found no consistent differences between the IUCN categories when compared deforestation rates and changes in the human footprint index [27]. Similarly, a second global study found that although deforestation rates tended to be lower in stricter PAs, this pattern was not Box 1. The six categories of PAs as officially defined by the IUCN based primarily on their management objective. Reproduced courtesy of IUCN from Dudley N. (Editor) (2008) [14].

Ia-Strict Nature Reserve
Areas set aside to protect biodiversity and also possibly geological/geomorphological features, where human visitation, use, and impacts are strictly controlled and limited to ensure the protection of the conservation values.

Ib-Wilderness Area
Large unmodified or slightly modified areas, retaining their natural character and influence, without permanent or significant human habitation, are protected and managed so as to preserve their natural condition.

II-National Park
Large natural/near-natural areas protecting major ecological processes, along with characteristic species and ecosystems, which also provide environmentally and culturally compatible spiritual, scientific, educational, recreational, and visitor opportunities. III-Natural Monument or Feature Areas set up to protect a specific natural monument, which can be a landform, seamount, submarine cavern, geological feature such as a cave, or even a living feature such as an ancient grove. IV-Habitat/ Species Management Area Set up to protect particular species or habitats with management reflecting this priority. Many but not all such areas will need regular, active interventions to meet the requirements of specific species or to maintain habitats.

V-Protected Landscape/ Seascape
Areas where the interaction of people and nature over time has produced an area of distinct character with significant ecological, biological, cultural, and scenic value, and where safeguarding the integrity of this interaction is vital to protecting and sustaining the associated values. VI-PA with sustainable use of natural resources Generally large areas, mostly in a natural condition, where a proportion of the total area is under sustainable natural resource management and where low-level nonindustrial use of natural resources compatible with nature conservation is seen as one of the main aims. consistent across all of the world's subregions [33]. Another global study, which examined the extent of croplands within PAs, found that there were more croplands in categories VI and V (followed by category IV) compared to categories I-III [10]. No study has yet analyzed the levels of human settlements [34] within PAs. The prevalent but untested assumption is that less strict PAs have higher levels of human settlements because of their management objectives (box 1). Although not all human presence and activities are incompatible with conservation [35], higher levels of human settlements could potentially translate into higher levels of human pressure on biodiversity [30,[36][37][38]. Consequently, PAs with more human settlements will necessitate conservation strategies that acknowledge and engage the local communities to successfully protect biodiversity [24,25,[39][40][41]. Our study's first objective was to assess whether the levels of human settlements within the world's PAs vary according to the categories, and particularly whether they tend to be higher in less strict PAs.
Our second but equally important objective was to identify the factors that may influence the levels of human settlements [34] within the world's PAs. Understanding which characteristics are associated most with higher levels of human settlements within PAs is essential for effective conservation. In addition to the IUCN categories [13], researchers have proposed a range of other possible drivers of human pressure within PAs [13]. However, it is not clear to what extent those factors influence human settlements and whether their influence varies across geographic regions. To address this second knowledge gap we evaluated the relationship between the following nine variables and the levels of human settlements [34] within the world's PAs: (a) elevation and (b) slope; previous studies have suggested that PAs at higher elevations and steeper slopes may have lower levels of human pressure, since areas with those characteristics tend to be less suitable for other human uses [42]; (c) bioclimatic zone; likewise, PAs located in areas with more extreme bioclimatic conditions [43] may have lower levels of human settlements because those conditions make them less attractive for human habitation [44]; (d) accessibility; PAs that are less remote, and closer to urban centers, are more accessible and therefore may have higher levels of human settlements, irrespective of the other characteristics [38]; (e) croplands; agriculture represents a major human economic activity, even within many of the world's PAs [10,12], and therefore it may serve as an important driver of human settlements; (f) country; human settlements' levels within PAs may vary according to the country in which they are located [9], reflecting many of the socioeconomic characteristics often suggested in the literature to influence the levels of human pressure within PAs, such as governance, management effectiveness, and availability of resources for conservation; (g) year of establishment; it is probable that older PAs have lower human settlements' levels since they have enjoyed longer periods of protection [10]; (h) area; previous studies have shown that smaller-sized PAs tend to have higher levels of human pressure [8]; (i) IUCN category; as we mention above, it is often argued that the levels of human pressure within PAs vary according to their management objective [17,19].

Data collection
We first downloaded the World Settlement Footprint (WSF) developed and made available by Marconcini et al [34]. WSF is a detailed map of the human settlements across the globe for the year 2015 and at a resolution of 10 m [34]. WSF represents mainly 'buildings' and 'building lots' and, to a lesser extent, 'roads and other paved surfaces' [34]. Although there are several other global maps of human settlements available [34], WSF is one of the most up-to-date, and is based on extremely high-resolution satellite information (10 m; figure 2). This is an important advantage for the type of analyses we performed in this study because it allowed us to use an exceptionally high number of PAs (i.e. 81 100-137 523 PAs compared to 12 315 PAs used in [9], 15 281 PAs in [33], and 19 486 PAs in [27]). To promote the use of WSF for global analyses, Marconcini et al [34] made available five resampled versions of the dataset, ranging from 100 m to 10 km (table S1). For the purposes of our analyses, we used the 100 m version, which essentially represents the percentage of human settlements within each 100 m × 100 m pixel (i.e. 0.01 km 2 ) across the globe.
To calculate each PA's proportion covered by human settlements, we used the World Database on Protected Areas (WDPAs), available online at protectedplanet.net (September 2020 version). Following the best practices guidelines, we removed all PAs not already established (i.e. those with 'proposed' status) and all UNESCO Man and Biosphere Reserves; these are essentially large areas in which buffer and transition zones are usually unprotected (http://protectedplanet.net/c/calculating-protectedarea-coverage). Additionally, we removed any PAs established after 2015 to avoid any temporal mismatches between the datasets used, which could result in inaccuracies. A small percentage of the PAs in WDPA (∼10%) does not have an establishment year; following previous studies, we randomly assigned to those PAs an establishment year based on the year that the rest of the PAs in the country were established [8,27]. To evaluate how this random assignment may have influenced our results, we re-ran our analyses using only the PAs that had already an establishment year. Both analyses produced nearly identical results, and therefore we report only the former, which is based on larger sample size. Lastly, we removed also all areas smaller than 1 km 2 to ensure that the WSF values calculated for each PA represented the human settlements within the PAs' boundaries (rather than their surrounding areas). A less conservative approach would have been to remove all areas smaller than 0.1 km 2 , which is still ten times larger than the resolution of WSF. Since both thresholds produced similar results, we only present below the analyses based on PAs ⩾ 1 km 2 (n = 81 100). The analyses based on PAs ⩾ 0.1 km 2 (n = 137 523) is presented in the supplementary materials (table S4 (available online at stacks.iop.org/ERL/16/064085/ mmedia)).
The final step of the data collection process involved extracting for each PA the following nine factors hypothesized to influence the human settlements' levels within PAs (a more detailed description of each factor is provided in table S1): (a) mean elevation and (b) mean slope, measured in meters and degrees respectively; (c) bioclimatic zone, identified using the global map of Metzger et al, who classified earth into 18 bioclimatic zones (table S1) based on a series of temperature and rainfall data, at a resolution of 1 km [43]; (d) accessibility, measured as the number of minutes required to travel from the nearest urban center to the PA [45]. We calculated accessibility using the 'global map of travel time to the nearest city' developed by Weiss et al available for the year 2015 and at a resolution of 1 km [45]. Weiss et al developed their global map using spatial data on transportation infrastructure (e.g. roads and railways) and other important landscape features known to affect human movements, such as topographical characteristics, land-cover types, and national borders. They defined as cities all areas with at least 1500 inhabitants per km 2 or areas with a population center with more than 50 000 people [45]; (e) extent of croplands, measured as the percentage of the PA covered by croplands and calculated using the Global Land Cover map (v 3.0) for the year 2015 available at a resolution of 100 m [46]; (f) country in which the PA is located and (g) year of establishment, as listed in the WDPA and calculated using the method described above; (h) area of each PA, measured in km 2 using the spatial boundaries of the PAs available in WDPA (table S1); (i) official IUCN Category, as reported by each country and listed in the WDPA.

Data analysis
To assess the extent to which the above nine factors are related to the levels of human settlements [34] within the PAs, we used the random forests (RF) technique, an ensemble machine learning method, which is based on decision trees [47,48]. Decision trees are used often to identify the factors that best predict a variable of interest-in this case the human settlements' levels. They are built by dividing incrementally a bootstrapped subset of the original dataset using a randomly selected predictor at each increment and evaluating its performance (table S2). Being an ensemble method, RF models combine the results of many such trees to create more accurate predictions [49]. We developed a separate RF model for each of the world's 14 subregions defined by the United Nations (figures 3 and 4(A)). The subregions are analogous to those used by the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services in their subregional assessments on biodiversity and ecosystem services [50]. Running a different model for each subregion was necessary in order to take into account regional variations in PAs [8,11] and to be able to decipher important regional patterns. Had we ran only one model, with all the PAs included, the results would have mostly reflected the patterns of the few subregions in the Global North that have a disproportionately large number of PAs (e.g. Northern America, Northern Europe, and Western Europe; table 1).
The RF method is used commonly to model global patterns [48,51,52] and has been shown to perform better than linear regression models [49] and to produce accurate results [53]. A key advantage of this particular method is that it does not make any assumptions regarding the relationships between the variables [54]. Moreover, RF is known to handle well large datasets with multiple predictor variables [54]. However, RF is an aspatial method [54], meaning that it does not take into account the spatial arrangement of the variables analyzed [48]. Consequently, and depending on the nature of the analysis performed, RF can sometimes produce spatially biased results [54]. A suggested approach to overcome this issue is to add the spatial coordinates of the units analyzed in the RF models [48,54]. The Moran's I test can be then used to evaluate whether the results are spatially biased by measuring the levels of spatial autocorrelation in the residuals of the models [48]. Moran's I values close to −1 or +1 indicate strong spatial autocorrelation (negative and positive respectively), while Moran's I values close to zero suggest no autocorrelation. Our preliminary analyses, following the above-described approach, showed that all 14 RF models produced near-zero Moran's I values (table 1); however, five of those values were statistically significant (table 1), mostly due to the large sample sizes used in the analyses, which provide the test with enough statistical power [55] to detect even minor levels of spatial autocorrelation (table 1).
To ensure that our results were robust and not affected by any issues associated with spatial autocorrelation, we used a second method to analyze each subregion's data, called Geographical Random Forests (GRFs). The method, which was developed by Georganos et al is essentially a spatial extension of the RF method. The idea behind GRF is similar to the Geographic Weighted Regression [48]; for each location in the dataset (i.e. each PA in our case), GRF computes a separate RF model using a userspecified number of nearby observations (table S5). Consequently, the resulting models are calibrated locally rather than across the whole subregion, producing less spatially biased results [48]. Since both methods (RF and GRF) produced similar results, leading to the same conclusions regarding our two research objectives, for brevity we only present here the results of the RF method. The results of the GRF method are presented in the supplementary materials (table S5; figures S9-S12).
To assess the relative importance of each variable included in the models, we used the two most common metrics for this purpose: (a) the percent increase in mean square error (%IncMSE) and (b) the increase in node purity (IncNodePurity) [56]. In both cases, the higher the value, the larger the variable's influence on the levels of human settlements. For each RF model, we also noted the percent of variance explained, which is essentially a measure of how well the model fits the data. Lastly, since the WSF valuesused as the response variable-represent a proportion (i.e. the proportion of each PA covered by human settlements), we applied the logit-transformation to ensure that the RF models' predictions were bounded between 0 and 1 [57].

Results
The mean percentage of human settlements within the PAs varied considerably between the 14 subregions examined (figure 3; table S3). Eastern Asia, Western Europe, and Southern Europe had on average the highest percentage of PAs' area covered by human settlements (7.5%, 2.1%, and 2.1% respectively; table S3), while Australia and New Zealand, and Northern America had the lowest (0.1% and 0.3% respectively).
The percentage of variance explained by the models also varied between the subregions (table 1). The highest percentages were in Eastern Asia, Latin American and the Caribbean, and Sub-Saharan Africa (80.6%, 65.7%, and 60.5% respectively), while the lowest were in Northern Africa and Pacific Islands (30% and 10.7%, respectively). Importantly, for at least eight of the fourteen subregions, our models explained more than 50% of the variation in the levels of human settlements within PAs (table 1).
There was also variation in the relative importance of the nine explanatory variables used in the analyses (figures 4(B) and S1-S4). That said, there were some consistent and important patterns. The first key pattern was that the differences in the levels of human settlements between the IUCN categories were minor in all 14 subregions (figures 4(B) and S5-S8). This was reflected by the low %IncMSE and IncNodePurity values of the specific variable in all the models, suggesting that the IUCN categories were a poor predictor of the human settlements' levels (figures 4(B) and S1-S4). The only exception was Table 1. Results of the RF models showing the number of PAs used in the analyses (n), the mean square error (MSE), the percentage of variance explained by each model (%var), the levels of spatial autocorrelation in the residuals (Moran's I value), and whether that autocorrelation was statistically significant (p-value). Significance levels: * p < 0.05, * * p < 0.01, * * * p < 0.001. Figure 5. Examples of the differences in the human settlements' levels between the IUCN categories in four of the subregions examined. Each data point represents a PA and each blue horizontal line represents the predicted value for that particular category according to the RF models. As it can be seen below, the levels of human settlements within PAs varied widely but the differences between the categories were, for the most part, minor. The results for the rest of the ten regions, which are nearly identical to those shown above can be found in the supplementary materials (figures S5-S8). XX = PAs with no IUCN category.
Western Europe in which the IUCN categories had the second-largest %IncMSE value and the fourth-largest IncNodePurity value (figure 6; table S6); this was mostly due to the somewhat higher levels of human settlements within PAs in category V compared to the rest of the categories ( figure 5). In the rest of the subregions, the categories' differences were largely insubstantial (figures 5 and S1-S4). Importantly, in Figure 6. Examples of the relative importance of the variables used in each model, as measured using %IncMSE and IncNodePurity. In most of the regions, accessibility was the most important predictor of the human settlements' levels within PAs.
The results for the rest of the ten regions can be found in the supplementary materials (figures S1-S4).
no subregion did we observe a pattern of increased human settlements from categories I to VI (figures S5-S8). The second key pattern, which was also consistent across all subregions, was that accessibility was associated strongly with the human settlements' levels within PAs (figures S1-S4). Specifically, in all but one of the regions, accessibility ranked the highest, both in terms of %IncMSE and IncNode-Purity (figures S1-S4). In fact, in many regions, such as in Southern Asia, Latin American and the Caribbean, and Sub-Saharan Africa, the relative importance of accessibility was considerably larger than any of the other variables included in the model (table S6; figures S1-S4). This essentially means that most of the variance explained in those models (table 1) was due to accessibility ( figure 4(B)).
Overall, the rest of the variables included in the analyses were less influential (table S6), and their importance varied across the subregions. For example, croplands were only important in Western Europe and Eastern Asia and somewhat important in Southern Europe, Northern Europe, and Southern Asia (table S6; figures S1-S4). Similarly, elevation and slope were only important in a few of the regions, such as in Eastern Asia and Southern Europe (figures S1-S4).

Discussion and conclusions
Our study is one of the first to analyze the patterns and potential drivers of human settlements within the world's PAs. The spatial scale of our analyses and its large sample size (n = 81 100-137 523) suggest that our results can be generalized across most of the world's PAs. Therefore, our findings could be used by policymakers and practitioners worldwide to improve the management of PAs. Moreover, as the parties to the Convention of Biodiversity Diversity prepare for their upcoming COP-15 meeting in October 2021-at which they will finalize the post-2020 global biodiversity targets, including those concerning the expansion of the global network of PAs [5]we hope that our findings can be used by governments to inform their plans and policies regarding the placement of future PAs.
A key finding stemming from our analyses is that the human settlements' levels within the world's PAs vary little between the IUCN Categories (figures S5-S8). Although it is reasonable for policymakers, conservation practitioners, and researchers to assume that strict PAs will have lower levels of human settlements, our results suggest that this is likely not the case and that the differences between the categories are minor in most parts of the world (figures S5-S8). Our findings are in agreement with the recent findings of Elleason et al, who showed that the changes in human footprint index [58] and forest loss across the IUCN categories were, for the most part, small and not statistically significant [27]. Multiple reasons could explain these patterns. First, it is important to note that the IUCN classification system is voluntary with no mechanism in place to verify that the PAs meet the conditions of their assigned category [16]. Therefore, there could be a considerable mismatch between the IUCN categories reported by each country and the realities on the ground. For example, in a recent analysis it was shown that several of Madagascar's PAs do not meet the conditions of their assigned category [59]. Second, it is important to note that there are numerous types of PAs worldwide, which do not always meet the exact criteria set by the IUCN [15,27]. As a result, there are likely to be considerable differences even between PAs that belong to the same category [27].
A second key finding of our analyses is that accessibility was one of the most predictors of the levels of human settlements within PAs. This pattern was consistent across all of the 14 subregions examined (figures S1-S4) and suggests that PAs nearer urban centers have higher human settlements' levels irrespective of their other characteristics (such as management objective, establishment year, and topographical features). Consequently, larger attention may need to be paid to those PAs. Managing successfully PAs with higher human settlements' levels will require conservation strategies that engage the local communities [26,60,61] and acknowledge the rights and the needs of the peoples in those areas [23,26,62].
Some of the other variables that may have influenced the levels of human settlements within PAs, at least in some of the subregions, include croplands, elevation, and slope. Croplands were associated with increased levels of human settlements, especially in Western Europe and Eastern Asia, and to a lesser extent in Southern Europe, Northern Europe, and Southern Asia (figures S1-S4). This suggests that agriculture could be an important factor influencing the levels of human presence within PAs, especially in those subregions. It must be noted, however, that agriculture could also be the result of higher levels of human settlements. Therefore, the relationship between these two variables needs to be explored further using a more appropriate dataset and study design. PAs located at lower elevations and on less steep slopes could also potentially experience higher levels of human settlements in at least some of the subregions (figures S1-S4).
The importance of the rest of the variables was consistently low (figure 4(B); table S6) suggesting no strong relationship with the levels of human settlements within PAs. To some extent, however, those patterns may have been influenced also by the uncertainty associated with measuring and capturing those relationships and processes at the global level. A somewhat surprising finding was that the country in which the PAs were located was rarely an important predictor of their human settlements' levels. We had expected that the country's identity would have had a larger influence because it captures many of the socioeconomic characteristics that are often mentioned in the literature as being important determinants of human pressure within PAs [9,13]. Examples of those characteristics include governance effectiveness and availability of resources for managing PAs, such as personnel, funding, etc [13]. In retrospect, country was not as important as originally hypothesized, most likely because our analyses was conducted at the subregional level (figure 4(A)), and at that level, many of the countries have very similar socioeconomic characteristics.
Having said that, it is worth noting that several recent studies, which have looked at the importance of some of those socioeconomic characteristics, have found the associations to be weaker than expected. For instance, governance effectiveness was not important in explaining the levels of human pressure within PAs at the national level [12], as measured using the human footprint index [58]; many of the countries that scored high on measures of national governance, such as several of the Western European countries, had also some of the highest levels of human pressure within PAs and vice versa [12]. Likewise, the association between conservation investments at the national level and levels of human pressure within PAs was also weak [12]. Analogous findings were found at smaller scales. For instance, Schleicher et al, who focused on 43 PAs in the Peruvian Amazon, measured the relationship between various management indicators (such personnel, funding, and infrastructure) and forest degradation and deforestation and found that the 'associations were at best weak' [63]. Similarly, a recent global study did not find any associations between management effectiveness, measured using Management Effectiveness Tracking Tool, and the levels of human pressure within 407 PAs worldwide [9].
Yet, as the authors of the above-mentioned studies explain the absence of a strong association does not necessarily mean that there is no relationship between management effectiveness and the levels of human pressure within PAs [12]. Quantifying and measuring management quality and effectiveness is challenging, especially at the scales analyzed in the above studies, and even more so here. Moreover, the unexplained variance in our models (table 1) suggests that other factors besides the nine tested are likely to be important-perhaps some of those factors include variables associated with management. We hope that as more data on management, including community participation, become available in the future [9], researchers will be able to analyze and understand better its influence on the levels of human settlements and pressure within PAs [63]. Nevertheless, despite our analyses not providing many insights about the importance of management, our main findings and conclusions, particularly regarding the IUCN categories and the importance of accessibility, remain valid and significant, considering especially: (a) the global coverage of our analyses, (b) the large number of PAs included, and (c) the large amount of variation explained for most of the 14 subregions (table 1).
We should clarify here that our analyses concerns only human settlements, as measured specifically using the WSF [34], and does not necessarily reflect biodiversity outcomes. Many would rightfully argue that human settlements, and human activities in general, are not necessarily incompatible with conservation. In fact, there are cases in which certain human activities, such as low-intensity farming, may be essential for the preservation of specific habitats and species [64,65]. Although we agree with this view, we believe that the information presented in this study has useful conservation implications, especially considering the current global rates of population growth and urbanization [38]. PAs with high and increasing levels of human settlements will require conservation strategies that engage the local communities to protect biodiversity successfully. Our findings could help governments and policymakers worldwide to focus their often-limited resources on the PAs that need them the most, and to plan more effectively the expansion of their network of PAs within the Post-2020 Global Biodiversity Framework [5].
To be clear, however, we do not advocate that governments take actions to reduce human settlements within PAs [66]. Those familiar with the vast literature on PAs and their social impacts [23,24,62] will know that there has been a long struggle to acknowledge local peoples' rights to continue their livelihoods within PAs [67]. It is now widely accepted that PAs must be managed equitably [24,68] and must not result in environmental injustices [69,70]. The current IUCN classification system is in a sense a response to these ideas [21,22,68] and it is for that reason that several of the categories (box 1) allow by design a range of sustainable nonindustrial activities [16,71]. Interestingly, though, our results suggest no clear relationship between the IUCN categories and the human settlements' levels. Instead, other factors, such as accessibility, appear to be better predictors of those levels. Policymakers and conservation practitioners will need to focus on those and will need to be ready to acknowledge and engage the local communities within PAs [26,60,61].

Data availability statement
All data used in the study are freely available online.
No new data were created or analyzed in this study.