Protected areas (PAs) are essential to prevent further biodiversity loss yet their effectiveness varies largely with governance and external threats. Although methodological advances have permitted assessments of PA effectiveness in mitigating deforestation, we still lack similar studies for the impact of PAs on wildlife populations. Here we demonstrate the application ofuse an innovative combination of matching methods and hurdle-mixed models with a large-scale and long-term dataset of unprecedented coverage for Finland’s large carnivore species. We show that the national PA network , at the national level, PAs does not support higher densities than non-protected habitat for 3 of the 4 species investigated. For the brown bear, PAs appear to have lower densities than non-protected areas. For some species, PA effects interact with region or time, i.e. wolverine densities decreased inside PAs over the study period and lynx densities increased inside eastern PAs. Although we show that matching approaches could and should be applied to wildlife population data, Wwe support their application of matching methods in combination of additional analytical frameworks for deeper understanding of conservation impacts on wildlife populations. These methodological advances are crucial for improving PA targets and extremely timely for preparing ambitious PA targets a post-2020 global framework for biodiversity.

Wildlife population time series data

Track count data collected within the Finnish Wildlife Triangle Scheme were used as indices of lynx, wolf, wolverine and bear densities. The scheme is a long term, large-scale survey which provides yearly estimates of the distribution and relative abundance of game species. The methods are described in detail by Linden et al.55. A wildlife unit is a permanent line transect route of 12 km (4 x 3 km). Unit locations represent different forest habitats in proportion to local occurrence. The monitoring operation is carried out twice a year. The minimum number of monitoring experts is three for the late-summer count and one for the winter count.

Each winter, 800–900 units distributed throughout Finland are surveyed by local hunters who count the number of fresh snow tracks crossing the transect. The count is either done 1–2 days after a snow fall or 24 h after a pre-check where all old tracks are marked. For each species, the data are compiled by the Finnish Game and Fisheries Research Institute as the number of crossings per 24 h per 10 km. Thus, the unit used in statistical modelling was the average number of crossings per 24 h per 10 km for lynx wolf and wolverine, per triangle and for each year. As brown bears hibernate during winter counts, density indices for this species were collected during late summer counts and also focused on counting tracks crossing the transect.

From the initial 2171 sampling units, a total of 1805 units were included in subsequent analyses (we discarded transects in mixed agricultural-forest landscapes), 1195 were located in non-protected sites while 610 were located inside protected areas.

Explanatory variable selection and extraction

The set of matching covariates extracted for each wildlife unit represented ecological characteristics considered most likely to be important determinants of large carnivore occurrence in European ecosystems based on empirical observation56. These covariates can be regarded as two groups of possible influences: biophysical context (e.g. latitude) and human impacts (e.g. distance to closest settlements). Several variables suggested by the literature to be important in determining PA effectiveness, e.g. PA management, were available for a restricted amount of sites and therefore were not included in the matching analysis. Description of the matching covariates and the rationale for their inclusion are found in the Supplementary Table S8.

Spatial data were analyzed using the R packages ”sp”, “raster”, “rgdal” and “rgeos”57. PA boundaries were calculated using spatial information from the World Database of Protected Areas58. A sampling unit was considered as protected if the 1km circular buffer centered on the unit’s centroid intersected a protected area polygon. Other explanatory variables corresponding to each wildlife unit were extracted for the unit’s centroid –for distance variables or gridded variables with a spatial resolution larger than the unit’s dimensions– or over a circular buffers centered on the unit’s centroid and passing through the three vertices –for gridded variables with a spatial resolution smaller than the unit’s dimensions (Supplementary Table S8).

Estimating protected area effectiveness with matching methods

Matching methods

We used the Matchit package59 in R environment, which fits a logistic generalized linear model where the treatment assignment (land protection) is the response variable and the matching variables are the predictors.

We chose one-to-one, nearest-neighbor covariate matching with replacement using a generalized version of the Mahalanobis distance metric. We used a caliper of 0.2 standard deviations of the propensity scores as our minimum matching criterion. To assess the quality of the matches, we checked the resulting covariate balance testing the normalized differences in covariate means and their distribution. The normalized difference in means is the difference in the average covariate value divided by its standard deviation60. We tested for differences in the distribution using eQQ plots that graph the covariate values in the same quantile of the treated (protected sites) against those in the control (non-protected sites), allowing us to observe if characteristics are distributed similarly across both treatment and control groups (see Supplementary Figures, S1a and S1b).

We were able to match 100% of the original sample to controls that suited the criteria. Therefore, our unit-to-unit matching yielded a final dataset of 1220 sampling units, i.e. 610 protected wildlife units that were matched to 610 similar unprotected units across Finland. These matched pairs of PA and non-PA sampling units covering, a period of ~30 years, were included in subsequent analyses. These sampling units were not homogeneously distributed among the four Game Management clusters with Central Finland being the most represented cluster and Lapland the least (see Supplementary Table S2).

Measuring the absolute effect of PAs

Previous studies have quantified the effect of PAs in reducing threats in different ways (Absolute PA effect, Relative PA effect, Pooled relative effect17). However, due to limitations imposed by the structure of wildlife time series and the natural low densities of large carnivores (zero inflation of all time-series), an effectiveness metric based on the ‘absolute PA effect’ was the most relevant approach for this study.

The absolute PA effect is the difference between densities in a unit located inside a PA and its matched control unit located in a non-protected area. Therefore a positive value means that sampling units located inside PA show higher large carnivore densities than its control unprotected units. We calculated this metric at the national PA network level, for each pair of ‘protected-non protected’ units. We computed the median absolute effect of the PA network and its associated 95% confidence interval for each species of large carnivore, globally and at the regional level. We performed iterative random sampling (1000 iterations) to control for differences in the number of repeated observations (number of years) among pairs of matched units. This was done by selecting only one value of absolute PA effect for each pair of matched units before estimating the median PA effect across units. To calculate the absolute PA effect per region, we pooled the 15 Game Management Areas covering the whole country in four main clusters: South= ’Satakunta’ + ‘Etelä-Häme’ + ‘Kaakkois-Suomi’ + ‘Uusimaa’ + ‘Varsinais-Suomi’ + ‘Pohjois-Häme’; Central= ‘Etelä-Savo’ + ‘Rannikko-Pohjanmaa’ + ‘Keski-Suomi’ + ‘Pohjois-Savo’ + ‘Pohjois-Karjala’ + ‘Pohjanmaa’; North= ‘Oulu’ + ‘Kainuu’; Lapland = ‘Lapland’. We compared the annual absolute PA effect per year with zero using one sample t-tests.

Assessing Estimating PA effectiveness with two-part mixed effects model for semi-continuous data.

We extracted the matched dataset obtained through the matching process described previously (1220 unpaired wildlife units, see Supplementary Table S7) to test for the effect of land protection status on large carnivore densities. We implemented two-part zero-inflated mixed effects models in the novel GLMMadaptive package that uses adaptive Gaussian quadrature61. This approach allowed us to account for the data structure (zero-inflated, right-skewed continuous data, GLMMzi) and assess the relationship between density indices of the four large carnivore species and the set of explanatory variables described in Table S8. All models were implemented in R using the packages GLMMadaptive and parallel. To account for potential problems of pseudoreplication, unit identity number was kept consistently as a random effect in all models for each of the 4 large carnivore species. Zero-inflated structures were added on all the fixed effects included in the models.

We reduced the full list of variables based on co-linearity and biological relevance to produce a set of 8 variables (e.g. collinearity between distance to closest roads and distance to closest settlements was too high and we chose to remove distance to closest roads from all the models). These 8 covariates included the 6 matching covariates described earlier (percentage of forest cover, terrain ruggedness, distance to closest settlements, human population density, latitude and longitude), to which we added year and protection status of the wildlife unit (protected or not, coded 0/1). Three interaction terms between fixed effects were also added (between PA and Latitude, PA and Longitude and PA and years) to assess spatial and temporal variations in PA effectiveness.

A set of 12 models was built for each species. All 12 models included the 6 matching covariates and protection status at the unit scale (PA), therefore the difference in model structure resided in the addition of year, the three interactive terms and their different possible combinations (structure of the 12 models is described in the Supplementary Table S9). All first-order model fits were ranked using the Akaike Information Criterion, the best model having the lowest AIC values from the set of 12 models built for each species.

The fixed effects estimates in mixed models with nonlinear link functions have an interpretation conditional on the random effects. To take this into account, we extracted the marginalized coefficients and their standard errors from the two part mixed models following the approach described by Hedeker et al.62 and implemented by Rizopoulos61 in the GLMMadaptive package.

Goodness-of-fit of the models was assessed using residual diagnostics following the procedures described in the DHARMa package. All statistical analyses were performed using the software R 3. 5. 157.

In order to test if PA size could affect our results, we built four additional models following the same procedure highlighted above. Data were restricted to protected wildlife units, density of the four species of carnivores was the response variables and covariates included: percentage of forest cover, human population density, distance to the closest settlements, terrain ruggedness, latitude, longitude, year and PA size. We extracted the marginalized coefficients and checked the residuals of the different models as highlighted above.