Temperature niche composition change inside and outside protected areas under climate warming

Conservation of biodiversity relies heavily on protected areas but their role and effectiveness under a warming climate is still debated. We estimated the climate‐driven changes in the temperature niche compositions of bird communities inside and outside protected areas in southern Canada. We hypothesized that communities inside protected areas include a higher proportion of cold‐dwelling species than communities outside protected areas. We also hypothesized that communities shift to warm‐dwelling species more slowly inside protected areas than outside. To study community changes, we used large‐scale and long‐term (1997–2019) data from the Breeding Bird Survey of Canada. To describe the temperature niche compositions of bird communities, we calculated the community temperature index (CTI) annually for each community inside and outside protected areas. Generally, warm‐dwelling species dominated communities with high CTI values. We modeled temporal changes in CTI as a function of protection status with linear mixed‐effect models. We also determined which species contributed most to the temporal changes in CTI with a jackknife approach. As anticipated, CTI was lower inside protected areas than outside. However, contrary to our expectation, CTI increased faster over time inside than outside protected areas and warm‐dwelling species contributed most to CTI change inside protected areas. These results highlight the ubiquitous impacts of climate warming. Currently, protected areas can aid cold‐dwelling species by providing habitat, but as the climate warms, the communities’ temperature compositions inside protected areas quickly begin to resemble those outside protected areas, suggesting that protected areas delay the impacts of climate warming on cold‐dwelling species.


INTRODUCTION
An increasing number of species are under threat from climate change and many populations are at risk of decreasing in abundance.Consequently, some species have become endangered (Amano et al., 2020;Bateman et al., 2020;Spooner et al., 2018;Tayleur et al., 2016).Therefore, predicting how the range, distribution, abundance, and phenology of species may change in response to a warming climate is a timely and important step in planning actions to mitigate the negative impacts of climate warming.In response to the changing climate, common generalist species or species with common functional traits tend to expand their ranges, which homogenizes community compositions (Davey et al., 2013;Fourcade et al., 2021).The distributions of many species are shifting toward higher elevations (Auer & King, 2014;Freeman et al., 2018) and latitudes (Chen et al., 2011;Fossheim et al., 2015) as a result of climate warming.However, range shifts are not constant across space, such that the leading range edges move faster than the trailing range edges (Massimino et al., 2015).Moreover, changes in species ranges tend to lag behind shifts in temperature (Devictor et al., 2008;Gaget et al., 2021;Santangeli & Lehikoinen, 2017;Savage & Vellend, 2015).Other global change drivers may also act in synergy with climate warming.For example, anthropogenic land-use changes can reduce species colonization in higher latitude and elevation habitats and aggravate declines of cold-dwelling species (Lehikoinen et al., 2019;Oliver et al., 2017;Virkkala et al., 2014).
Many species are more dependent on protected areas (PAs) than on the typically more human-modified landscapes surrounding them (Coetzee et al., 2014;Gillingham et al., 2015;Virkkala et al., 2014).Earlier studies show that PAs can facilitate colonization by providing habitats (Thomas et al., 2012).
Moreover, habitat depletion, species' extinctions, and population declines can be lower inside PAs (Geldmann et al., 2013;Peach et al., 2019).In general, PAs play an important role in conservation and in maintaining biodiversity, whereas intensive land use continues to decrease the quality and quantity of unprotected areas (Geldmann et al., 2019;Joppa et al., 2016).PAs can mitigate the negative effects of climate change on population trends (Gillingham et al., 2015;Lehikoinen et al., 2019).Indeed, species' abundances remain higher inside PAs than outside (Gillingham et al., 2015;Pavón-Jordán et al., 2020).For instance, the abundance of northern species tends to be higher inside than outside PAs (Lehikoinen et al., 2019).In general, the negative effects of a warming climate are less pronounced in protected natural areas compared with unprotected areas surrounding them, where anthropogenic pressures are often rampant (Häkkilä et al., 2018;Northrup & Rivers, 2019).PAs reduce the loss of natural cover nearby when land surrounding a PA is conserved as a private nature reserve, indicating that the PA border is not a strict ecological border (Ament & Cumming 2016 ).Moreover, species do not heed PA borders; thus, community-level spillover may occur from PAs to the matrix (Boesing et al., 2021).
Despite the positive effects of PAs for biodiversity, their ability to buffer biodiversity loss in the face of climate change is contested (Monzón et al., 2011).To stay in their current climatic conditions, species may need to shift their distributions.However, most PAs have been established to conserve species and habitats in specific locations, but due to climate-driven range shifts, the targeted species and habitats may not benefit from those static PAs in the future.The boreal network of PAs are biased toward higher elevations (Sanderson et al., 2015), where climate warming alters conditions more rapidly (Pepin et al., 2015).Therefore, the contemporary network of PAs may not be effective and able to prevent the negative impacts of climate change on species (Holsinger et al., 2019;Johnston et al., 2013;Wu et al., 2018).However, analyses of the impact of PAs in mitigating climate change effects on ecological communities' at large spatial scales are few.In large-scale studies, PAs can be considered a network, rather than individual sites.Moreover, studying PA networks at a large scale allows the assessment of PA mitigation of climate change effects on ecological communities.
The ability of PAs to mitigate climate-driven changes in species' abundances and occurrences can be studied using the community temperature index (CTI), which quantifies the average temperature niche of a community.The CTI has been used to study climate change responses of several taxonomic groups (Devictor et al., 2008;Fourcade et al., 2019;Savage & Vellend, 2015).The CTI measures the balance between cold-and warmdwelling species (cold-dwelling species prefer areas with colder temperatures compared with warm-dwelling species) in a community (Devictor et al., 2008).Thereby, high values of CTI reflect a community dominated by warm-dwelling species and vice versa.Combined with spatiotemporal occurrence data of species, CTI can be used to explore whether the thermal signature of ecological communities is shifting in response to climate warming (Devictor et al., 2012a).
The CTI combines species abundances and traits and may thus be driven not only by climatic forces, but also by other drivers of population change, such as land-use change (Clavero et al., 2011).However, several studies have documented increases in CTI (Curley et al., 2022;Lehikoinen et al., 2021) and mismatches in the rates of change of CTI and temperature change (Bertrand et al., 2011;Santangeli & Lehikoinen, 2017).Nevertheless, it is often unclear whether a CTI change is due to an increase in the abundance of warm-dwelling species or a decrease in the abundance of cold-dwelling species.For example, in fragmented landscapes, habitats without strong anthropogenic pressure favored warmdwelling butterfly species (Fourcade et al., 2021), but a strong anthropogenic pressure amplifies the declines of cold-dwelling bird species (Oliver et al., 2017).The CTI change can also be a result of simultaneous change in abundances of cold-and warm-dwelling species (Tayleur et al., 2016).We aimed to fill this knowledge gap by disentangling the contributions of warmand cold-dwelling species to the overall change in temperature niches in bird communities.
We studied community-level responses of birds to warming climate inside and outside of PAs at large spatial scales over time.For this, we used roadside survey data from the North American Breeding Bird Survey (BBS) scheme for Canada from 1997 to 2019.We explored whether PAs influence the speed of climatedriven changes in communities' temperature niches, quantified as CTI.More specifically, we asked does overall CTI of bird communities differ inside and outside PAs, does the temporal trend in CTI vary inside and outside PAs, and which species contribute most to the difference in the temporal change of CTI between areas inside and outside PAs?
For the first question, our hypothesis was communities inside PAs have a higher proportion of cold-dwelling species than areas outside PAs because high-latitude PAs located in boreal ecosystems aim to conserve habitat for cold-dwelling species, whereas unprotected areas may be more suitable for warm-dwelling generalists (Lehikoinen et al., 2019;Santangeli & Lehikoinen, 2017).We expected that PAs would more strongly mitigate the effects of climate change for cold-dwelling species than for warm-dwelling species because the land cover inside PAs changes less relative to unprotected lands, where land-use change and other anthropogenic pressures are increasing (Leroux & Kerr, 2013).We did not study the effects of habitat type; rather, we assumed, based on earlier studies, that PAs with certain characteristics, such as dense forest cover, can maintain high abundances of habitat specialist and cold-dwelling species (Barnagaud et al., 2013;Gaüzère et al., 2016;Lehikoinen et al., 2019).
For the second question, we expected that the positive temporal changes in CTI have occurred inside and outside PAs (Curley et al., 2022;Gaget et al., 2021;Santangeli et al., 2017) because the ambient average temperatures increased during the study period (Devictor et al., 2008;Lindström et al., 2013).However, we considered two competing hypotheses about differences in CTI change between PAs and unprotected lands.First, CTI change is slower in PAs.Based on previous studies, areas inside PAs can have less pronounced negative effects on the abundance of cold-dwelling species (Lehikoinen et al., 2019), which would lead to slower CTI change in PAs.This pattern could be because PAs contain natural areas (Fourcade et al., 2021), such as old-growth forest (Häkkilä et al., 2018), and suitable microclimate for cold-dwelling species (Frey et al., 2016).Second, CTI change is faster inside PAs because warm-dwelling species utilize PAs when colonizing new areas, thus increasing the CTI in PAs (Thomas et al., 2012;Fourcade et al., 2021;Gaget et al., 2021).
For the third question, in terms of species' contributions to the temporal changes in CTI inside and outside PAs, we considered two competing hypotheses.First, cold-dwelling species have a greater effect on the CTI inside than outside PAs, following recent findings (Gaüzère et al., 2016;Lehikoinen et al., 2019).This could cause reduced CTI changes compared with unprotected areas.Second, warm-dwelling species speed up the increase in CTI because abundances of warm-dwelling species increased in North America during the study period following recent climate warming (Curley et al., 2022;Princé & Zuckerberg, 2015).

Bird survey data
We used data from the North American BBS to assess the CTI changes in bird communities over time.The North American BBS is a standardized annual roadside point count survey that provides long-term and large-scale population information for hundreds of North American bird species (USGS Patuxent Wildlife, 2022).Altogether, over 4800 routes have been established using a spatially stratified and random design across the United States, Canada, and northern Mexico, although we only used data from Canada.For each route surveyed, a skilled observer conducts 50 3-min point counts evenly spaced along a 40 km section of road (point counts approximately 800 m apart).The survey is conducted during the peak of the breeding season (typically in June), and the observer records the species and number of individual birds heard from any distance or seen (except dependent young) within 400 m of each roadside point count location.We excluded all non-native species from the analyses because their distributions are less likely than those of native species to be at climatic equilibrium and because the species' climatic preferences are more likely to be unreliable when calculated for an unstable, human-driven, distribution (Dyer et al., 2017).We did not exclude any native species because the aim was to study the overall communities at the study sites.
We used the BBS data from Canada for 1997-2019 for which geographic coordinates of point-count locations were available.We assumed that point-count locations remained constant throughout the study period.The availability of fine-scale geographic coordinates of each survey point allowed us to link the survey point locations with the Canadian PA network.To allow assessment of temporal changes in the average temperature niches of bird communities, we used routes that had at least 2 years of survey data during the study period (see detailed route selection in Appendix S1).This resulted in a final data set of 821 routes, largely distributed across southern Canada (Figure 1).
On average, each route was surveyed 12 times (SE 0.24).The final data included 383 species and ∼6.1 million observed individuals.We used the geographic coordinates of each survey point to classify the protection status of the land at each point-count location (outside or inside PA).We delineated an area of 400-m radius around each survey point with GIS tools in ArcGIS 10.8 (ESRI, 2020).If the area overlapped a PA, the survey point was classified as protected, following, for example, Terraube et al. (2020) and Santangeli et al. (2022).We used a 400-m radius because it ensured that the surrounding landscape that affects observed species was considered and because the landscape just outside of PAs may also be influenced by the PA through spillover effects, such as a reduced threat to natural cover (Ament & Cumming, 2016;Shen et al., 2022).
Surveys are conducted along roads that are not necessarily protected even though the area alongside the road is protected.Often routes pass through protected and unprotected areas (286 routes).We split the point locations in each route into two groups: those inside PAs and those outside PAs.To reflect their differing locations along the 40-km route, we calculated group-wise average coordinates by averaging the coordinates of the points falling inside or outside PAs.If the route fell completely inside (nine routes) or outside (526 routes) PAs, we calculated the average coordinates for the route based on all points.Thereby, each study unit (hereafter, average site) was defined by the route and the inside or outside PA status of the points.For each average site, we summed up all observed individuals by species.We aggregated the data to avoid random variation of single 3-min point counts.
We obtained the Canadian PA network data from the World of Database on Protected Areas (IUCN, 2021), including all PAs in the International Union for the Conservation of Nature categories I-VI (Dudley, 2008).In these categories, the management restrictions range from strictly PAs to areas with sustainable use of natural resources, but nature conservation is one of the main aims in all of them.The establishment year of the PAs varied, but we assumed that sites that were designated later in the study period were of good environmental quality.Although PAs cover parts of all terrestrial Canadian ecozones (Environment & Climate Change Canada, 2016) and are distributed broadly, the network is unevenly distributed in space.The larger and more remote PAs tend to be in northern Canada, where land use is less intensive (Venter et al., 2016) than in southern Canada, where the PAs are smaller but more numerous.
We used climate data to examine the spatial variation in average temperature niches of bird communities inside and outside PAs.To account for the general climatic conditions in the study area, we acquired monthly averaged temperatures from weather data for 1970−2000 across Canada at 30-arc seconds resolution (WordClim database; Fick, 2017).From the weather data, we calculated the breeding season (March-August) mean temperatures for each average site.If the coordinates of the site did not overlap the temperature data layer, we obtained the mean temperature value within a 3-km area around the average site.

Community temperature index
To measure how the average temperature niche of bird communities has changed in space and time, we calculated CTI for each bird community for each year in each average site.The CTI is a community-weighted mean of species-specific temperature niches and derived from the species' climatic preferences (measured using species temperature index [STI]) (Devictor et al., 2008(Devictor et al., , 2012)).The STI represents the average temperature experienced by a species within its geographical range for a given season, reflecting its association with temperature.We used STI values calculated in Lehikoinen et al. (2021).They calculated the breeding season STIs based on species' breeding ranges according to Handbook of the Birds of the World and BirdLife International Illustrated Checklist of the Birds of the World (del Hoyo & Collar, 2016) and matched the ranges of species with the mean temperatures of the typical breeding months (March-August) from 1950 to 2000 (Hijmans et al., 2005).They restricted the breeding ranges to North America (including Canada, the United States, Mexico, and Greenland).The STI is not an absolute measure of the species temperature niche because it does not distinguish the broadness of the niche.Rather, STI should be considered a relative measure (i.e., an index) of a species' temperature affinity (Devictor et al., 2012b(Devictor et al., , 2012a)).
We summed all the individuals of the same species for each average site.After that, we calculated the CTI for each bird community in each average site in each year with breeding bird season STI values and species' counts.A high CTI value represents a community dominated by warm-dwelling species, whereas a low CTI value represents a community dominated by cold-dwelling species.

Long-term changes in CTI and contributions of different species
To study how the average temperature niches of bird communities have changed over time inside and outside PAs, we used a linear mixed-effects model with a Gaussian distribution.We considered the average site-specific CTI the response variable and the continuous year (standardized from 1997 to 2019), protection status (protected, unprotected), mean temperature of the average site (continuous), and an interaction between year and protection status as predictor variables.Initially, we also added the longitude and latitude to the model, but these variables were excluded because they were strongly positively correlated with each other and with mean temperature.The mean temperature and the protection status did not correlate with each other (Appendix S1).
We included route and grid cell identity as random effects to account for the spatially nested structure of the data.More specifically, we established a grid of 5 • × 5 • covering all of Canada, and we assigned each average site to a grid cell.By doing this, we were better able to attribute the observed variation in CTI to the variables of interest, rather than the spatial proximity of the protected and unprotected average sites.We weighted the models with the log-transformed number of points with which the CTI was calculated.That is, we gave more weight in the analysis to those average sites with a larger number of survey points and thus lower uncertainty in the CTI estimates.In spatial correlograms, we found no signs of spatial autocorrelation in the residuals of the final model (Appendix S1) (Zuur et al., 2009).Moreover, we tested the sensitivity of the months that were used to calculate the STI and found no effect of the month selection on the results (Appendix S1).
We used a jackknife procedure (Crowley, 1992) to investigate which species contributed most to the modeled CTI changes.In particular, our interest was in identifying whether the warmor cold-dwelling species were driving the interaction effect of year and protection status on CTI.For that, we removed species from the data set one by one to iteratively recalculate the CTI, after which we reran the linear mixed effect model.For each species, we calculated the difference between the interaction coefficient from the global model and the jackknife model.A positive difference between model coefficients indicates that a species contributed by strengthening the relationship between year and protection (i.e., the interaction effect), whereas a negative difference indicates that a species contributed by weakening the interaction effect.
We transformed the obtained contributions into absolute values and log-transformed them to meet the requirements of normality.We also calculated the difference between each species' STI value and the average STI value of all species considered.A high relative STI value indicates that the species is warm-dwelling, whereas a low relative STI value indicates the species is cold-dwelling.Using the transformed contributions and the relative STI-values, we fitted a linear model with a quadratic effect to study whether the warm-or colddwelling species contributed most to the observed interaction effect coefficient.A high species' contribution suggests that the species is particularly influential in defining the difference in temporal CTI change inside and outside PAs.We did not correct the models with phylogenetic relationships among species because we were interested in the absolute effect that each species had on the interaction effect,  a The interaction coefficients are obtained from the model of community temperature index as a function of interaction between year and protection and mean temperature.b Relative STI describes the average temperature across the species' breeding range in relation to the average STI across all species in the data set.
regardless of the phylogenetic relationships within studied communities.
All analyses were carried out with R (R Core Team, 2019).For the linear mixed effect models, we used the functions lmer (Bates et al., 2015) and lmerTest (Kuznetsova, 2017).For the linear model, we used the function lm (Bates et al., 2015).Spatial autocorrelation was tested with package ncf (Bjornstad, 2022).

RESULTS
Patterns of average CTI and CTI change differed inside and outside PAs.The CTI was significantly lower inside than outside PAs across our study area during the study period (Table 1& Figure 2).As expected, areas with higher mean temperatures also had higher CTI values.We also found that CTI increased over the 22-year study period inside and outside PAs; the increase was faster increase inside PAs than outside (Table 1).At the start of the time series in 1997, CTI was substantially lower inside PAs.However, the difference in CTI between areas inside and outside PAs had almost disappeared by 2019 (Figure 2), due to faster increases in CTI inside PAs.
Species varied in their contributions to the interaction effect coefficient of CTI change over time between PAs and outside areas.There was high variability in the species' contributions that were suggestively associated with the species' relative STI values.For example, the savannah sparrow (Passerculus sandwichensis) (STI 5.3) had the highest contribution to the interaction effect coefficient and the king rail (Rallus elegans) (STI 18.9) had the lowest contribution to the interaction effect coefficient.There was little evidence that those species with STI values close to the average STI across species (i.e., relative STI ≈ 0) were contributing more than species with strongly positive or negative relative STI.Most importantly, there was a significant, but weak, linear trend (β = 0.02, SE = 0.007, 95% CI 0.006-0.034,p = 0.027) that species with a high relative STI value (i.e., warm-dwelling species) contributed more strongly to the coefficient of the interaction between year and PA status (higher increases inside PAs).

DISCUSSION
Our results show that the average temperature niches in bird communities have changed at different rates inside and outside PAs in southern Canada over 22 years.Although the average temperature niche was consistently warmer outside PAs, the change in the average temperature niche was faster inside PAs.
Lower CTI values inside PAs than outside PAs have also been in Europe (Gaüzère et al., 2016;Santangeli et al., 2017).For example, in northern Europe, PAs are often located in boreal old-growth forests (Auvinen et al., 2010), which have a cooler microclimate during the avian breeding season (Frey et al., 2016), which makes these PAs more suitable for cold-dwelling, boreal biome specialists (Santangeli et al., 2017).This means that PAs can act as a climatic refuge for northern species as the climate warms.However, climatic conditions are not unequivocally the primary factor determining geographic distributions (Rich & Currie, 2018) or the main driver of distribution shifts (Platts et al., 2019).The quality of habitat in the boreal forest is also associated with higher species density (Häkkilä et al., 2018) and lower CTI (Lehikoinen et al., 2019).The considerable range of CTI changes in our results suggests that habitat quality or other mechanisms may have influenced our results.For example, higher proportions of old-growth forest can preserve greater species richness (Häkkilä et al., 2018), and different timber harvesting methods may change bird abundances (Cadieux et al., 2020).
In accordance with our hypothesis and inconsistent with some earlier studies (Santangeli et al., 2017), the temperature niches of bird communities changed faster inside PAs than outside PAs.That is, the temperature niches of bird communities inside PAs largely resemble those of the bird communities outside PAs.However, the level of uncertainty was high and further studies examining the habitat type and quality variation inside and outside PAs could potentially explain some of the uncertainty.That said, the faster change inside PAs may be the result of the PAs ability to facilitate species' range expansions (Thomas et al., 2012;Peach et al., 2019) rather than faster declines of cold-dwelling species (Lehikoinen et al., 2019).Indeed, our results show that warm-dwelling species contribute to the overall change in CTI in the bird communities more than cold-dwelling species, which may suggest that they utilize PAs more effectively as stepping stones to expand to new areas and thus contribute to the overall change in CTI in the bird communities than cold-dwelling species.Similar patterns have been observed elsewhere, especially in fragmented landscapes, where warm-dwelling species have exhibited higher colonization probability relative to cold-dwelling species (Fourcade et al., 2021).PAs are often placed at higher elevation (Sanderson et al., 2015), and the climate warming rate seems to increase with elevation (Pepin et al., 2015).Faster warming could drive faster CTI change inside PAs than outside PAs because areas outside PAs tend to be located at lower elevations.
The warm-dwelling species may have a greater influence on the decreasing differences in bird communities' temperature niches inside and outside PAs based on our results, but the high uncertainty in the results prevents definitive conclusions.The effect of warm-dwelling species is supported by PAs ability to promote the colonization of warm-dwelling bird species during the breeding season (Thomas et al., 2012).In general, colonization-extinction dynamics can be inconsistent across species' ranges, such that the leading edges of the ranges tend to expand more than the trailing edges of the ranges contract (Massimino et al., 2015;Oliver et al., 2017).Therefore, in the northern hemisphere in particular, finding habitat in a warming climate may be easier for warm-dwelling species than cold-dwelling species.It is also possible that warm-dwelling FIGURE 3 The relationship species' absolute contributions to the interaction effect between year and site protection status and relative species temperature index (STI) (point, a species removed in the jackknife analysis; red, positive effect on the interaction effect coefficient; blue, negative effect on the interaction effect coefficient).Relative STI is the difference between species-specific temperature index and the average temperature index across all species in the data set (i.e., North American Breeding Bird Survey, Canada from 1997 to 2019).Separate linear models with quadratic terms are shown for the positively and negatively contributing species for illustrative purposes.For the actual fitted model, see Table 2.
species have reached a certain abundance level outside PAs and, therefore, are expanding to PAs as a spillover effect.A similar spillover effect is driving generalist species from managed to natural areas (Frost et al., 2015;Häkkilä et al., 2017).Consistent with the spillover effect, earlier studies show that the CTI change in eastern North America is driven by warmdwelling species (Curley et al., 2022) and that the colonizations of warm-dwelling species drive the increase in the number of bird species overwintering in urban areas in North America (Princé & Zuckerberg, 2015).In parallel, cold-dwelling species are declining in northern European breeding bird communities (Tayleur et al., 2016).However, the fact that the association between the species' contribution to the interaction coefficient and relative STI values is hump-shaped rather than u-shaped (Figure 3) suggests that CTI changes are largely driven by abundance changes of species whose temperature niches do not largely differ from the average temperature niches of all species and not by abundance changes of species with extreme temperature niches.
Across our study area, even the most remote PAs are under some level of human influence, such as global climate change.Some of the largest and most remote PAs of Canada are in the northern parts of the country, where larger changes in climate have occurred than in the southern parts (Zhang et al., 2019).Moreover, roadside PAs across the country are becoming more isolated because of intensifying land-use practices (Andrew et al., 2011;Leroux & Kerr, 2013).In Canada, the roadside PAs are small and tend to have high primary production, but they are surrounded by intense land use, such as clearcuts (Andrew et al., 2011).High-intensity land use can be related to rapid declines in cold-and warm-dwelling bird species (Cadieux et al., 2020;Oliver et al., 2017), and highly fragmented PAs may not be able to preserve species richness (Durán et al., 2016).
To understand the context-dependency of how PAs maintain biodiversity under ever-intensifying climate warming, we encourage future researchers to account for finer-scale details in the environment surrounding the PAs, as well as the properties of the PAs themselves (Bowgen et al., 2022;Clavero et al., 2011).First, we considered a binary protection status of the survey points, but further nuances in that regard could provide information on potential thresholds in PA size or habitat heterogeneity that are required to buffer the bird communities against climate warming effects.We assumed that even the smallest amount of PA inside the 400-m radius around the survey point would have an effect on the bird community and that all the PAs maintained equivalently higher bird populations.Although this is a reasonable assumption given that birds are mobile and thus likely to use the PAs near the locations they were observed, we recommend future research on a finer scale to understand how much the size of the nearby PAs affects bird abundances.
Second, we considered recently established and old PAs, but their effects on bird communities may differ.The amount of protected land increased during our study period (Environment & Climate Change Canada, 2016).Therefore, it may be relevant to assess how the age of the PA affects the changes in the bird communities and whether the recently established PAs provide similar buffering against climate warming impacts as the old ones.
Third, based on earlier studies, we assumed that habitat types and quality of the habitats inside PAs are higher than outside PAs for cold-dwelling species (Häkkilä et al., 2018;Lehikoinen et al., 2019), but we did not account for fine-scale land use.Because there was notable uncertainty in CTI changes outside and inside PAs, we recommend future researchers to explore how local land use, habitat type, and quality or distance to the edge of the PA may influence changes in community composition (Sirami et al., 2017).The last limitation is that we did not control for observer experience, which could cause a skewness in the observed species such that large-bodied and more easily detectable birds could be overrepresented.That said, detectability is a bigger issue with unstructured citizen science data (Callaghan et al., 2021), but we used systematically collected data.
Our results support the role of PAs in mitigating the impacts of climate warming on birds, but raise a warning that the extent of this mitigation may be limited if temperatures continue rising.We suggest that climate warming should have a more prominent role in PA management strategies.For now, the challenges brought by the warming climate have been taken into account, but the strategies should be updated to enhance adaptations to intensifying climate change (Barr et al., 2021).Alarmingly, Canada and the United States have reduced their governmental support for development and management of PAs (Watson et al., 2014), and an increasing number of threats, such as introduced invasive species (Rose & Hermanutz, 2004) and loss of old-growth forest (Määttänen et al., 2022), add pressure on the PAs in the boreal zone in general.The fact that communities inside PAs are under the same threats as the ones outside PAs, but the communities' responses are different, highlights the need for additional research at finer scales.To bridge the gap between research and PA management, researchers should explore the fine-scale context-dependencies of climate change effects on biodiversity in PAs.

FIGURE 1
FIGURE 1The averaged coordinates of the grouped breeding bird survey routes in southern Canada based on the point's protection status (i.e., the average sites in Canada from 1997 to 2019) (gray open circles, averaged unprotected sites; black crosses, averaged protected sites).

FIGURE 2
FIGURE 2 Change in community temperature index southern Canada from 1997 to 2019 in areas inside and outside protected areas (shading, 95% CI).

TABLE 1
Summary results of the linear mixed-effect model explaining long-term change in community temperature index in Canada from 1997 to 2019.a a The category protected was used as the reference category of protection status.

TABLE 2
Summary of the quadratic regression model of bird species' contributions to the interaction effect coefficient as a function of relative species temperature index (STI).a