Assessing the effectiveness of a national protected area network for carnivore conservation

The study addresses whether Finland’s national protected area (PA) network is ecologically effective at conserving large carnivores by supporting higher population densities than comparable non-protected areas. While PAs are widely considered key conservation tools, empirical evidence on their effectiveness for wildlife populations remains limited, especially when accounting for non-random PA placement and confounders. Building on advances in counterfactual evaluation (matching) used for deforestation, the authors investigate PA impacts on the densities of Eurasian lynx, gray wolf, wolverine, and brown bear across Finland over ~30 years. They hypothesize that effective PAs should yield higher carnivore densities, but predict variation among species (with wolverine expected to benefit most), across regions (potentially lower effectiveness in Southern Finland with small PAs and in Lapland due to reindeer-related conflicts), and over time (possible declines in effectiveness with increasing conflicts). The work aims to inform PA policy and management by providing a rigorous, wildlife-focused effectiveness assessment.

Prior research has emphasized PA efficiency (e.g., range representation) over ecological effectiveness, with growing calls to evaluate outcomes for biodiversity. Counterfactual approaches (matching) have shown that PAs reduce deforestation, though impacts are often overestimated without accounting for location biases. For wildlife, some studies reported more positive trends inside PAs or stable trends in tropical PAs, and contributions of PA networks to megafauna habitat, but these often lacked counterfactual controls, risking overestimation. Large carnivores globally face habitat loss and persecution, with European populations recovering in parts of the continent; in Fennoscandia, legal hunting occurs under quotas but illegal killing is a major concern, sometimes even within PAs. In Finland, adaptive management and increased prey have aided recoveries, though wolf growth has been constrained by poaching. These contexts underscore the need for rigorous, counterfactual-based evaluations of PA impacts on carnivore populations.

Study system and data: The Finnish Wildlife Triangle Scheme (FWTS) provides long-term, nationwide carnivore abundance indices from ~12 km permanent line transects (triangles of 4 × 3 km). From 2171 initial units, 1805 were analyzed (610 in PAs; 1195 outside), excluding mixed agricultural–forest transects. Winter counts (800–900 units annually) recorded fresh snow tracks for lynx, wolf, and wolverine, producing indices as number of crossings per 24 h per 10 km per unit-year. Brown bear indices were collected during late summer due to hibernation.

Defining protection and covariates: A unit was considered protected if a 1 km circular buffer around the centroid intersected a PA polygon (WDPA). Explanatory variables captured biophysical context and human impacts: percentage forest cover, terrain ruggedness, distance to closest settlements, human population density, latitude, and longitude. Spatial extractions used R packages sp, raster, rgdal, and rgeos. Variables such as PA management were unavailable at sufficient coverage and thus excluded.

Matching to estimate absolute PA effect: Using MatchIt in R, the authors performed one-to-one nearest-neighbor matching with replacement based on a generalized Mahalanobis distance and a 0.2 SD caliper on propensity scores. Balance was assessed via normalized mean differences and QQ plots. They matched 610 protected units to 610 similar unprotected controls (1220 units), spanning ~30 years and distributed across four regional clusters (South, Central, North, Lapland). The effectiveness metric was the absolute PA effect: for each matched pair and year, the difference in density between the PA unit and its matched non-PA unit. They computed the median absolute PA effect and 95% CIs per species at national and regional scales, using iterative random sampling (1000 iterations) to select one observation per pair to account for unequal time-series lengths. Annual effects were compared to zero via one-sample t-tests. Regional analyses pooled 15 game management areas into four clusters.

Hurdle mixed-effects models (two-part zero-inflated GLMMs): To corroborate and probe spatio-temporal heterogeneity, they modeled densities using GLMMadaptive with adaptive Gaussian quadrature, accounting for zero-inflation and right-skew. The response was the density index per unit-year, with unit ID as a random effect. Fixed effects included the six matching covariates plus year and PA status (0/1). Interactions tested spatio-temporal variation in PA effects: PA × latitude, PA × longitude, PA × year. Collinearity led to excluding distance to roads (retaining distance to settlements). For each species, 12 candidate models (all included the six covariates and PA; differing by inclusion/combinations of year and interactions) were ranked by AIC. Marginalized coefficients and SEs were extracted following Hedeker et al. and implemented in GLMMadaptive. Model diagnostics used DHARMa. Additional models restricted to protected units tested PA size effects, including the same covariates plus PA size.

Sample sizes: From matched data, species-specific paired samples included lynx (n=207 units), wolverine (n=75), wolf (n=67), and bear (n=209). The hurdle models used 1126 unpaired sampling units (post-matching dataset preparation) covering the national extent and study period (1989–2017).

• National-level matching (all years, all regions) showed no significant PA effect on lynx, wolverine, or wolf densities: lynx median absolute PA effect −0.077 (95% CI: −0.250, 0.091; n=207); wolverine 0.097 (95% CI: −0.250, 0.091; n=75); wolf 0.018 (95% CI: −0.031, 0.077; n=67). Brown bear densities were lower inside PAs: median absolute PA effect −0.311 (95% CI: −0.523, −0.081; n=209).
• Annual analyses found no year-specific differences between PAs and non-PAs for any species (all Bonferroni-adjusted p-values = 1 for years with adequate observations).
• Regional matching analyses (South, Central, North, Lapland) did not reveal significant spatial contrasts (all 95% CIs overlapped 0).
• Hurdle-mixed models corroborated no global PA effect on lynx (estimate 0.0041; z=0.624; p=0.533), wolverine (−0.0018; z=−0.659; p=0.509), and wolf (0.0008; z=0.330; p=0.742). For bears, hurdle models indicated higher densities inside PAs (estimate 0.0065; z=1.875; p=0.0607), contrasting with matching results.
• The bear discrepancy likely reflects dataset structural differences (fewer paired samples in matching) and scale: matching estimates pair-level differences, whereas mixed models estimate marginalized effects across all data; spatial heterogeneity may allow a few high-density protected sites to drive a positive global effect.
• Interactions revealed spatio-temporal heterogeneity: • Lynx: PA × longitude interaction significant (estimate −0.018; z=−5.567; p=0.010), with higher densities outside PAs in western Finland but higher inside PAs in eastern Finland. • Wolverine: PA × year interaction negative (estimate −0.005; z=−2.029; p=0.042), indicating stable densities outside PAs but declining densities inside PAs over ~30 years; effect disappears when excluding Lapland, suggesting Lappish sites drive the decline.
• PA size had no significant effect on densities for any species (bear: 0.0011; z=0.269; p=0.788; lynx: 0.0019; z=0.698; p=0.484; wolf: −0.0032; z=−0.949; p=0.342; wolverine: −0.0003; z=−0.104; p=0.917).

The analyses indicate that Finland’s PAs do not generally support higher densities of lynx, wolf, or wolverine compared to similar non-protected areas, with bears showing method-dependent signals. This suggests complex interactions among law enforcement, hunting/poaching pressures, and resource distribution. High ungulate abundance and productive anthropogenic landscapes outside PAs may offset potential benefits of reduced mortality or better habitat within PAs for carnivores reliant on ungulates, explaining limited PA effects for lynx, wolf, and wolverine. Bears’ omnivorous diet may reduce dependence on ungulate-rich non-PA landscapes, contributing to the different pattern.

Spatial and temporal heterogeneity matters: lynx benefit more from PAs in the east, potentially due to prey (mountain hare) availability, high habitat suitability, and proximity to Russian source populations. Conversely, wolverine densities declined inside PAs over time, particularly in Lapland, aligning with heightened human–carnivore conflict and potential increases in anthropogenic mortality (including poaching and licensed depredation control) within protected areas in reindeer husbandry regions. Seasonal reindeer use of old-growth resources within PAs might elevate conflict and mortality risk inside PAs in the north.

These findings emphasize that lack of overall density differences does not imply PA irrelevance but highlights governance, enforcement, conflict mitigation, and landscape context as key determinants of outcomes. Methodologically, combining matching with mixed models offers complementary insights—pair-level counterfactual contrasts and population-level marginalized effects—critical for interpreting PA impacts on mobile wildlife.

This study provides, to the authors’ knowledge, the first application of statistical matching to wildlife population data, combined with hurdle-mixed models, to evaluate PA effectiveness for large carnivores at national scale. Overall, Finland’s PAs did not yield higher densities for lynx, wolf, or wolverine relative to matched non-PAs, with method-dependent signals for bears and important spatial (east–west) and temporal (declining wolverines inside PAs in Lapland) heterogeneity. The work underscores the need for targeted management, improved co-management and enforcement in conflict-prone northern PAs, and consideration of prey dynamics outside PAs. Methodologically, it advocates using counterfactual matching alongside flexible mixed models to capture heterogeneity and strengthen inference. Future research should: incorporate direct measures of anthropogenic pressures and PA management quality; use fine-scale GPS data to assess seasonal habitat use; evaluate connectivity and transboundary dynamics; and generalize the framework to other taxa and regions to inform post-2020 PA targets and SDG 15.

• Equivalence of densities between PAs and non-PAs does not necessarily indicate ineffectiveness; unmeasured pressures (e.g., enforcement, poaching, legal culling) and benefits (habitat quality, prey) may counterbalance.
• Seasonal mismatch in data collection (bears surveyed in late summer, other carnivores in winter) could bias PA effect detection due to seasonal habitat selection and prey dynamics.
• Use of track counts as indices, while validated for these species, entails indirect density estimation and potential measurement error.
• High mobility of large carnivores and movement across PA boundaries complicate attributing effects to protection status.
• Limited availability of key covariates (e.g., PA management quality, enforcement) restricted inclusion in models, potentially confounding results.
• Matching reduces sample size to paired units and may miss broader patterns; mixed models rely on model specification and assume appropriate zero-inflation structure.
• Spatial heterogeneity may lead to contrasting conclusions across methods (e.g., brown bear), highlighting sensitivity to dataset structure and scale of inference.