Non-native Species Surrounding Protected Areas Influence the Community of Non-native Species Within Them

Protected areas (PAs) are a key element of global conservation strategies aiming to protect habitats and species from various threats such as non-natives species (NNS) with negative ecological impacts. Yet little is known about the mechanisms by which PAs are colonized by NNS, and more specifically the role of colonizing events from surrounding areas. Here, we compared terrestrial and freshwater non-native plants and animals recorded in Norwegian PAs and in 5-km belts around them, using the database of the Norwegian Biodiversity Information Centre Species Map Service. Our analysis included 1,602 NNS and 671 PAs. We found that NNS were recorded in only 23% of the PAs, despite the fact that 90% of the 5-km belts were colonized by at least one NNS. A Zero-inflated negative binomial regression model showed that the number of NNS in the 5-km belts was a strong explanatory variable of the NNS richness inside PAs. Other significant variables included the surface area of the PA, mean human population density in the PA, main type of habitat and accessibility of PAs. We also observed similarity in the species in and around the PAs, with, on average, two thirds of the NNS present in a specific PA also present in its 5-km belt. Furthermore, NNS were recorded in PAs on average 4.5 years after being recorded in the 0–5 km belts, suggesting a dynamic of rapid colonization from the belts to the PAs. Invasive NNS represented 12% of NNS in the belts but 40% in the PAs. This difference was related to the higher abundance of invasive NNS in the belts. Our results highlight the necessity of expanding the focus of NNS management in PAs beyond their boundaries, in particular to prevent incursions of NNS with high negative ecological impact.


INTRODUCTION
Protected areas (PAs) are key target elements of global biodiversity conservation strategies. In 2010, 150 governmental leaders committed, through the Aichi Target 11 in the Convention on Biological Diversity, to improve the status of biodiversity by setting 17% of the global terrestrial area under protection by 2020 (CBD, 2020). The main purposes of PAs are to maintain natural ecosystem functioning, prevent habitat degradation due to human activities (Rodrigues et al., 2004), conserve biodiversity (Worboys, 2015) and protect nature from various threats (Mathur et al., 2015) such as acting as natural filters against invasive non-native species (Foxcroft et al., 2011).
According to IUCN, a non-native species (NNS) is a species introduced outside its natural past or present distribution (IUCN, 2016). During the last two centuries the number of introduced NNS species has increased substantially worldwide with no sign of saturation (Seebens et al., 2017a). Their spread is a consequence of increased human mobility, and the expansion and globalization of trade between countries and continents (Nunes et al., 2015;Chapman et al., 2017;Seebens et al., 2018;Ward et al., 2020). Although the ecological impacts of most NNS are either negligible or unknown (Jarić and Cvijanović, 2012;Seebens et al., 2018;Blackburn et al., 2019), some non-natives are invasive, or potentially invasive: i.e., they have negative impacts on the recipient species and ecosystem (IUCN, 2016;Blackburn et al., 2019). Biological invasions are one of the leading causes of global biodiversity loss (Intergovernmental sciencepolicy platform on biodiversity ecosystem services, 2019) and are one of the principal drivers of recent species extinctions (Clavero and García-Berthou, 2005;Bellard et al., 2016;Blackburn et al., 2019).
Numerous guidelines and technical tools have been developed to assist in the management of invasive NNS in PAs (e.g., de Pooter et al., 2007;Monaco and Genovesi, 2014). These manuals generally advocate the early detection and eradication of all NNS, including those that have not been proven to be invasive, as an implementation of a precautionary approach that considers all NNS to be potentially invasive (McNeely et al., 2001;Monaco and Genovesi, 2014). Beyond the threat to biodiversity posed by invasive species, all NNS represent a human footprint on natural environments. NNS introduced by humans are considered undesirable in PAs, the purpose of which is to preserve nature in as pristine a state as possible (Hettinger, 2001). The presence of NNS also potentially contributes to increasing homogenization of native biological communities (McKinney and Lockwood, 1999;Lambdon et al., 2008;Kortz and Magurran, 2019).
Previous studies have shown that NNS richness patterns within PAs are linked to anthropogenic factors such as road networks and human population density inside PAs (Spear et al., 2013;Dimitrakopoulos et al., 2017;Gallardo et al., 2017;Moustakas et al., 2018). Other properties of the PAs, such as their surface area and protection status, also influence NNS richness (Gallardo et al., 2017;Liu et al., 2020). Furthermore, NNS presence in PAs is also driven by the properties of the surrounding areas, including human land use, human population density and road density (Foxcroft et al., 2011;Spear et al., 2013). These results suggest that, even if long-distance dispersal can be important for the expansion of NNS, especially in the early stages (Ramakrishnan et al., 2010), short-distance dispersal represents a significant contribution to their colonization dynamics. One of the few documented examples was published recently by Liu et al. (2020), based on the global alien distributions of 894 animal species: they found that 89-99% of PAs had an established population of at least one of these species within 10-100 km of their boundaries, but the majority of PAs were not colonized by any of them. Nevertheless, little is known about the influence of the NNS pool present in close proximity to PAs on the NNS communities within them (but for an example see Meiners and Pickett, 2013).
Here, we analyze the composition of terrestrial and freshwater non-native plants and animals present in Norwegian PAs, and in belts of 0-5, 5-10, and 10-20 km around them, to assess the extent to which the community of NNS in areas immediately surrounding PAs relates to NNS within PAs. We selected Norway due to the availability of an extensive database on NNS from the Norwegian Biodiversity Information Centre (NBIC) and the Global Biodiversity Information Facility Norway (GBIF Norway). We hypothesized that the presence of NNS in PAs should mainly be a result of colonization from surrounding areas. NNS in close proximity to PAs should thus influence the community of NNS present in PAs qualitatively, quantitatively and temporally. We expected to find: 1. A high proportion of NNS present in a PA are also present in its surroundings (qualitative similarity). The NNS in a PA will have taxonomic and ecological similarities to the pool of NNS in its surroundings. However, since invasive NNS are expected to have a higher colonization potential than non-invasive NNS, invasives should be present in higher proportions inside PAs than outside in comparison to other NNS. 2. The total number of NNS present inside a PA is a positive function of the richness of NNS in its surroundings (quantitative influence). In addition, the most abundant species in the surroundings of PAs are more likely to be present within the PAs. 3. NNS are recorded in the areas surrounding PAs earlier than inside the PAs (temporal sequence).

Data on Non-native Species
We downloaded NNS data from the NBIC Species Map Service (https://www.artskart.artsdatabanken.no, 10/04/2020. Data from: List supplementary material. Downloaded through the Species Map service). This database is provided by various contributors including research institutes, environmental agencies and NGOs. Biodiversity data from online databases are potentially biased, for example by accessibility of sites, lack of coverage of geographic and environmental variation that cover species distributions (Hortal et al., 2007), or by taxonomy, such as societal preferences in citizen science projects (Troudet et al., 2017). However, we consider the Norwegian database as one of the most robust that is available, as it contains an extensive collection and evaluation of NNS from a wide range of taxa and across the country (Sandvik et al., 2019;Tsiamis et al., 2019). We selected terrestrial and freshwater NNS records of the Kingdom "Animalia" and "Plantae" with an accuracy of <= 100 m. For this purpose, we retained only those species with the following habitat categories assigned by Norwegian Biodiversity Information Center (NBIC, https://www.biodiversity.no): terrestrial, limnic/terrestrial, limnic/marine habitats. We filtered for records from the year 1950 to the present. After selection, our NNS database included a total of 350,286 records of 1,602 species representing 21 different taxonomic classes. 14.9% of the NNS were animals and 85.1% were plants.
In order to consider potential ecological impacts caused by NNS, we used the ecological risk assessment conducted by the Norwegian Biodiversity Information Centre, in which each NNS is assigned to one of the following categories: -"Severe impact" (SE): NNS with actual or potential ecologically harmful impact and the potential to become established across large areas; -"High impact" (HI): NNS with either a moderate ability to spread but which cause at least a medium ecological effect, or have a minor ecological effect but have a high invasion potential; -"Potentially high impact" (PH): NNS with either a high ecological effect and low invasion potential or high invasion potential without known ecological effect; -"Low impact" (LO): NNS with no substantial impact upon Norwegian nature -"Not known impact" (NK): NNS with no known impact; -"Not risk assessed" (NR): NNS not yet risk assessed. NNS belonging to "Severe impact" and "High impact" categories are included in the Norwegian Black List 2012 of Alien Species. In total, 60% of the NNS included in the analysis were risk assessed, while 40% were not.

Data on Protected Areas
We extracted the shape files and information on the designation year and surface area of Norwegian PAs from the World Data Base on Protected Areas (WDPA, UNEP-WCMC and IUCN, 2019). The WDPA contains 3,143 registered Norwegian PAs of which 2,178 PAs are terrestrial and cover 54,749 km² of the land area of Norway (http://protectedplanet.net, accessed March 2020). We selected for analysis PAs with status "designated" and categorized as "terrestrial, " excluding "marine, " and "coastal" PAs. There are also PAs that are not assigned to any management categories (i.e., category marked as not assigned, not reported, not applicable); these PAs were excluded.
Protected areas of the WDPA are categorized in different International Union for Conservation of Nature (IUCN) categories, which range from category I (strictly protected or large, unmodified or slightly modified areas) to VI (protected areas with sustainable use of natural resources) whilst further PAs that are not assigned, not reported and not applicable (IUCN, https://www.iucn.org). For our study we selected PAs in category I, II (national parks) and IV (habitat/species management areas) with a surface area >= 1 km². Since our analysis investigated NNS in PAs and belt zones up to a distance of 20 km around the PAs (see below), we also excluded PAs whose belt zones crossed the political borders with Sweden, Finland and Russia.
Applying these filters resulted in 671 PAs in our analysis: 623 PAs of IUCN category I (average surface area = 7.8 km²), 18 PAs of IUCN category II (average surface area = 1,064.9 km²) and 30 PAs of IUCN category IV (average surface area = 142.3 km²). All PAs were designated between 1959 and 2017 and had areas ranging between 1 and 3,444.8 km² (average 42.2 km²). They covered 28,314.9 km², which is 49.5% of the total terrestrial protected area of Norway.

Belt Zones Around Protected Areas
We mapped belt zones of 0-5, 5-10, and 10-20 km circumjacent to PAs using QGIS (http://qgis.osgeo.org, version 3.4.2-Madeira) (Figure 1). All PAs and belt zones were entirely within Norway. Where the belt zone of a PA included part or all of another PA, the intersecting area was not excluded from the belt, such that belts should not be considered as indicators of the state of protection. Our analysis focused on the belt of 0-5 km (henceforth referred to as 5-km belt) to investigate whether the composition of NNS communities within PAs was influenced by NNS in the immediate vicinity of PAs. The surface area of this belt naturally varied with the size of the PA, with a range 99.1-2,218.2 km² (average 167.7 km²).

Statistical Analysis
All analyses were performed using the statistical analysis tool R [R Core Team (2019), https://www.R-project.org, version 4.2.0]. We extracted NNS records for all the PAs and their belts, deriving from this a list of NNS for each PA and its surrounding belt. To test qualitative similarity, we used tests across all PAs and belts (i.e., overall comparisons considering independently the records made in a PA and in 5-km belts). The temporal sequence analysis was based on data where the NNS was present in the PA and the associated belts using a pairwise comparison. To test for quantitative influence, we used a mixture of tests considering data in PAs and their associated belts as independent (tests on abundance) or as paired for the other analysis (NNS present in PAs and associated 5-km belts and modeling NNS richness).

Qualitative Similarity Taxonomy and Ecological Impact
We compared the proportions of taxonomic classes and ecological impact categories of NNS between the PAs and belts using Pearson's chi-squared tests.

Most Frequent NNS
To identify the most frequent NNS in the PAs, we selected NNS that were present in at least ten PAs. We compared them with the same number of NNS that were most frequent in the 5-km belts.

NNS Abundance
We defined the abundance of a species in a PA or a belt as the number of records of this species. The mean abundance of a species is therefore the number of records divided by the number of PAs or belts where it was present.
For NNS present in both the PAs and the 5-km belts, we tested if there was a correlation between their mean abundance in the PAs and the belt using a Spearman's rank correlation (rho).
To test whether the NNS present in the PAs are among the most abundant in the 5-km belts, we compared the mean abundance of NNS present in both the PAs and the belts with the mean abundance of NNS present only in the belts using a Wilcoxon Rank Sum Test.
To test whether NNS of the impact categories "Severe impact" and "High impact" (black listed NNS) were more abundant in the 5-km belts than NNS of less severe impact categories (nonblack listed NNS), we compared the mean abundance of these two groups of NNS in the 5-km belts using a Wilcoxon Rank Sum Test.

NNS Present in PAs and Associated 5-km Belts
To investigate the hypothesis that the presence of NNS in PAs is mainly a result of colonizing events from surrounding areas, we calculated, for each NNS present in the 5-km belts, the proportion of times it was present in both the belt and its associated PA, and the proportion of times it was present only in the 5-km belt but not in its associated PA. We then calculated the mean of these proportions for all the NNS present in the 5-km belts.
We applied the same approach for the NNS present in the PAs. We calculated the proportion of time they were present in both the PA and its associated 5-km belt, and the proportion of time they were present only in the PA but not in its associated 5-km belt. We then calculated the mean of these proportions for all the NNS present in the PAs.

Modeling NNS Richness
We selected five explanatory variables to model NNS richness in PAs, comprising two biotic variables (the most abundant land cover in the PAs and NNS richness in the 5-km belts); two anthropogenic variables (mean human population density of the region in which the PA is located and their accessibility), and PA surface area. Land cover was obtained from Copernicus Land Monitoring Service using information Label 1 (CLC, 2018). Label 1 information consists of five categories: "Artificial surfaces", "Agricultural areas", "Forest" and semi natural areas", "Wetlands" and "Water bodies". We extracted the land cover of each PA in QGIS and calculated the percentage of the most abundant land cover category in each PA.
NNS richness in each of the PAs and their associated 5-km belt zones was extracted from the NNS lists. The mean accessibility of PAs, calculated as the mean travel time from within PAs to the nearest city with a population >50,000 inhabitants, was extracted from Nelson (2008), a map integrating transportation networks and agglomeration index (a measure of urban concentration) and was downloaded from the European Commission (https://forobs. jrc.ec.europa.eu). Mean human population density was obtained from WorldPop (2018). The surface areas of the PAs were filtered from the World Database of Protected Areas (WDPA, https:// protectedplanet.net). We assessed the relationships between predictor variables using Spearman's Rank correlation. All predictors were retained for the analysis since they had little correlation among them (Supplementary Figure 1).
We applied a Zero-inflated Negative Binomial regression model with Poisson distribution (ZINB, Lawal, 2012) to test if the NNS richness in PAs is a result of the NNS richness in the associated 5-km belts, anthropogenic and PA properties. We assumed that if NNS have the opportunity to colonize PAs their richness inside is between 0 or higher and therefore is a count process (Gallardo et al., 2017). On the other hand, for PAs that were uncolonized by NNS, we assumed this was due to missing vectors, distance, or the PA not having suitable habitat (Gallardo et al., 2017). The only outcome in this case is zero. The ZINB model consists of two parts: The first part is the negative binomial regression model, which explains the relationship between conditional variance and conditional mean compared to the Poisson distribution model. The second part, the binary distribution model, captures the excess of zero values that exceed the predicted zeros by the negative binomial distribution. We used the package "pscl" to run the ZINB (Achim et al., 2008;Jackman, 2020).

Year of First Record in PAs and Associated Belts
To test whether NNS were recorded earlier in the surrounding belt than inside PAs, we extracted the years of first record for each NNS in the PAs and the 0-5, 5-10, and 10-20 km belt zones. We selected only NNS present in PAs. For each NNS we looked at the PA and the associated belts. We compared the years of first record in PAs and their three associated belts using a Kruskal-Wallis and Dunn's test for pairwise comparison with the "fdr" adjustment method.

RESULTS
We analyzed 671 Norwegian PAs, of which only 22.8% were colonized by any NNS. In contrast, at least one NNS was present in 89.5% of the 5-km belts. The total number of NNS records was 8,641 in PAs, and 156,736 in the 5-km belts, which represents 2.4 and 44.7%, respectively, of all the Norwegian NNS records included in the analysis. The remaining records were in the 5-10 km belts and 10-20 km belts or outside of them. The number of NNS was between 0 and 53 in PAs (mean = 0.87, SD = 3.68) and 0 and 440 (mean = 23.3, SD = 50.62) in 5-km belts. Of the 1,602 NNS in our analysis, 196 were present in the PAs, and 1,123 in the 5-km belts. All bar one of the 196 NNS present in the PAs (99.5%) were among the NNS present in the 5-km belts, the exception being the plant, Leucanthemum maximum. The number of NNS present in PAs also varied between IUCN categories (category I: mean ± SD = 0.68 ± 2.85; category II: 1.27 ± 2.35; category IV: 4.53 ± 11.03).

Qualitative Similarity
Taxonomy More than 75% of the NNS in both the PAs and the 5-km belts were plants, although the proportion of plants was lower in PAs than in 5-km belts (Figure 2A). Five plant classes were present in the PAs vs. eight in the 5-km belts ( Figure 2B). Eudicots (e.g., broadleaf trees Acer pseudoplatanus and Sambucus racemosa) represented the highest proportion in both the PAs and the 5km belts but the proportion was significantly lower in the PAs than in the 5-km belts ( Figure 2B). The inverse relationship was observed for Pinopsida (e.g., coniferous trees Picea stichensis and Abies alba), with a higher proportion of Pinopsida present in PAs. Non-native animal species represented 22% and 13% of the NNS in PAs and 5-km belts, respectively. Six animal classes were present in the PAs and 8 in the 5-km belts, of which Insecta (e.g., the beetles Acrotrichis insularis and Cartodere nodifer) showed the highest proportion in PAs and 5-km belts, followed by Aves (e.g., the Canada goose Branta canadensis and the Mandarin duck Aix galericulata), with a significant higher proportion of Aves in PAs ( Figure 2B).

Ecological Impact
NNS listed in the 2012 Norwegian Black List comprised ∼ 40% of the NNS in PAs, with 28.5% classified as species with "Severe impact" (SE) and 11.3% with "High impact" (HI) (Figure 2C). In contrast, 12% in the 5-km belts were listed in the Norwegian Black List, with 6.7% "SE" and 5.3% "HI, " this difference being significant (X² = 90, df = 1, p < 0.05). Fiftyseven percent of the NNS in the Black List that were present in the 5-km belt were also present in the PAs, compared to only around 12% of the non-listed NNS (X² = 165.19, df = 1, p < 0.001).

Most Frequent NNS
Nine NNS were present in at least 10 PAs. The most frequent being the Canada goose (Branta canadensis), which was present in 34/671 PAs (5%) (Figure 3). Of the top 9 NNS, six were plants and three were animals, one plant and one animal being aquatic. The three animals were chordates (B. canadensis, Neovison vison and Salvelius fontinalis). B. canadensis was also among the top 9 NNS present in the 5-km belts, but at a much higher proportion, colonizing 28% of them (Figure 3). In the 5-km belts, the most frequent NNS was the Garden lupin (Lupinus polyphyllus), which was present in 437/671 (65%) of the belts and was also among the top 9 NNS in PAs. Seven of the top 9 NNS present in PAs and seven of the top 9 NNS present in the 5-km belts were in the Norwegian Black List of Alien Species 2012.

NNS Abundance
The mean abundance of NNS present in PAs was significantly positively correlated with that of the 5-km belts (Spearman's rank correlation rho: S = 598,446, p < 0.001, rho = 0.52).
NNS present in both the 5-km belts and the PAs were significantly more abundant in the 5-km belts than the NNS present only in the 5-km belts (mean abundance in the belts: 8.51 and 2.13 records per NNS, respectively, Wilcoxon Rank Sum Test: W = 34,914, p < 0.001). In the 5-km belts, the abundance of black-listed NNS was significantly higher than the other NNS (mean abundance in belts: 9.48 and 2.38 records per NNS, respectively, Wilcoxon Rank Sum Test: W = 26,002, p < 0.001).

NNS Present in the PAs and Their Associated 5-Km Belts
Of the pool of NNS present in the associated 5-km belt of a PA, only 1% on average were also present inside the PA they surround. In contrast, on average 63% of NNS present in a PA were also present in its associated 5-km belt, while the remaining 37% of NNS present were only in the PAs and not in the associated 5-km belts.

Modeling NNS Richness
The results of the Zero-Inflated Negative Binomial regression show a significant positive relationship between NNS richness in PAs and richness in the associated 5-km belts (Table 1,  Figure 4). NNS richness in the 5-km belts was also the only significant variable in the zero part of the model (i.e., modeling PAs free of NNS), being lower when surrounding PAs with no recorded NNS.
The most abundant land cover in the majority of PAs was "Forest and semi natural areas" (525 PAs, 78.6%) followed by "Wetlands" (99 PAs, 14.8%) and "Waterbodies" (42 PAs, 6.3%) (Figure 4). Agricultural area was the most abundant land cover of only two PAs. The count part of the ZINB model shows that the number of NNS in PAs was highest where water bodies were the most abundant land cover ( Table 1). The number of NNS in PAs also significantly increased with increasing mean human population density in the PAs and the surface area of PAs, and was negatively correlated with travel time to large cities (Table 1,  Figure 4).

Year of First Records in PAs and Associated Belts
Overall, NNS were recorded later in PAs than in any of the three associated belts, with the difference on average being 4.5 years (0-5 km), 6 years (5-10 km), and 5.5 years (10-20 km) (Figure 5). The average years of first records in the three belts were not significantly different. Of the NNS present in both the PA and the associated 5-km belt, 59.4% were recorded earlier in the belt, 17.5% in the same year and 23.1% earlier in the PAs. This overall pattern of delayed records in the PAs was observed for 5 of the 11 taxonomic classes (Figure 6). Of the top 9 NNS in PAs, six were recorded significantly earlier in the PAs than in at least one of the belts (Kruskal-Wallis-Test: X² = 62.21, df = 3, p < 0.001) (Supplementary Figure 2).

DISCUSSION
Our study provides an extensive nationwide analysis of how the NNS community in the vicinity of PAs influences the NNS community inside PAs. Using data on 1,602 non-native terrestrial  and freshwater animals and plants of Norway, we showed that 77% of the PAs included in our analysis were free from any of them. This result is in accordance with a previous study on a global scale, which found more than 90% of PAs free from any of 894 non-native animals (Liu et al., 2020). The absence of NNS in PAs is often attributed to their remoteness, which keeps them far from areas where many NNS are introduced: the introduction of NNS is often associated with trading and transport activities between cities and countries (Banks et al., 2015;Nunes et al., 2015;Seebens et al., 2017b;Seebens, 2019) and NNS further spread by vectors such as roads, streams or intended and unintended human transportation (Leuven et al., 2009;Nunes et al., 2015;Brancatelli and Zalba, 2018;Ward et al., 2020). However, in our study, the low NNS richness in the 5-km belts surrounding PAs was the only variable explaining variation in the presence or absence of NNS in PAs (i.e., the zero part of the ZINB regression). This suggests that low colonization and propagule pressure in close proximity seems to be a better explanatory factor for the absence of NNS in PAs than their accessibility. For PAs occupied by NNS, their NNS richness was again significantly related to NNS richness in the surrounding 5-km belts, the accessibility of PAs having a lower effect (i.e., the count part of the ZINB regression). These results again support our hypothesis of a quantitative effect of the pool of NNS in areas close to PAs on the richness of NNS within the PAs. Nevertheless, three other factors also influenced the richness of NNS in PAs: their surface area, the human population density inside them and the main type of habitat they contain. PAs in which water bodies were the most abundant habitat had the highest NNS richness, highlighting lakes and rivers as corridors for the colonization of both limnic (Leuven et al., 2009) and terrestrial NNS (Malíková and Prach, 2010;Francis et al., 2019).
The temporal analysis carried out in our study revealed that NNS were recorded earlier in the immediate surroundings of PAs than within them. On average, NNS were recorded in the PAs 4.5 years after being recorded in the 0-5 km belts. We also measured a delay in the first records of NNS in the PAs for FIGURE 5 | Year of first records in protected areas and their associated belts 0-5, 5-10, and 10-20 km for the 196 non-native species recorded within protected areas. We selected only non-native species present in protected areas. For each non-native species, we compared protected areas and their associated belts. The numbers above the boxes represent the numbers of first records of the non-native species included in the analysis. Small letters (a, b) indicate elements that are significantly different from each other according to a Dunn's test for pairwise comparison following a significant Kruskal-Wallis Test (KW).
five of the eleven taxonomic classes of NNS and six of the nine most frequent NNS found in PAs. This spatio-temporal sequence of occurrence confirms that PAs are not prime locations for FIGURE 6 | Boxplot NNS species taxonomic classes present in PAs: The year of first record in PAs and within a distance of 0-5, 5-10, and 10-20 km from the PA (associated belts). The numbers above the boxes represent the numbers of first records of the non-native species. Small letters (a, b, c) indicate elements that are significantly different from each other according to Dunn's test with "fdr" adjustment, following a significant Kruskal-Wallis Test (KW). the introduction of NNS and further suggests the important role of colonizing events from within a few kilometers of their boundaries in the processes involved in the spread of NNS in PAs. These results also suggest an invasion debt, i.e., the time lag between the introduction of a non-native species into a region and its potentially negative ecological consequences (Rouget et al., 2016). For example, ornamental plants already introduced for horticultural purposes, but not yet naturalized (i.e., not yet established as persistent wild populations outside of cultivation), represent a risk of invasion in the future that could be exacerbated by climate change (Haeuser et al., 2018). Once naturalized, the time it takes for an invasive species to reach remote PAs, potentially containing many threatened native species, may be another element of this debt. Garden lupin (L. polyphyllus), for example, considered a severely impacting NNS, which was most common in the 0-5 km belts, but much less common in PAs, should require special consideration in their management. This pattern is also supported by the qualitative similarity that we observed within and around PAs, with, on average, two thirds of the NNS present in a specific PA also present in its associated 5-km belt. Nonetheless, previous studies have shown that successful colonization of new environments by NNS varies from species to species depending on environmental conditions and species characteristics (Sakai, 2001;Gallien and Carboni, 2017). Differences in environmental conditions inside and outside PAs, as already shown by Mas (2005), could explain differences in species frequency inside and around PAs. For instance, the Garden lupin (L. polyphyllus) the most frequent NNS in the 5-km belts, is an ornamental plant which is common in Norwegian gardens from where it escaped from cultivation (Fermstad, 2010). Another example of differences in environmental conditions are reflected by the fact that the proportions of NNS of two classes, Aves and Pinopsida, were higher inside PAs than in their surroundings. Ten of the 12 non-native birds were Anseriformes (ducks, geese and swans), thus dependent on aquatic habitats, such as the Canada goose (B. canadensis), which is the most frequent NNS in the PAs. The Canada goose utilizes open and grassy habitats and nearby lakes and other water bodies, feeding on aquatic plants and animals amongst other food (Jansson et al., 2008). Concerning the conifers (class of Pinopsida), suitable habitats comprise forest and semi natural areas, which was the most abundant land cover type in the majority of the PAs in our study. Non-native waterfowl and conifers may thus find more suitable habitat in PAs than around them, as many PAs, unlike belts, have probably been delineated to include high conservation value habitats such as water bodies and forests.
Our analyses show that invasive NNS (i.e., listed on the Norwegian, 2012 blacklist) are over-represented in Norwegian PAs compared to non-invasives. Invasive NNS accounted for 12% of the NNS in the 5-km belts but 40% in the PAs. Furthermore, 57% of invasive NNS present in 5-km belts are also present inside PAs. This high colonization success of invasive NNS in PAs may be explained by their high abundance outside PAs and by having characteristics that permit their fast colonization and spread. In the belt, an invasive NNS was, on average, four times as abundant as a non-invasive NNS (with species abundance measured as the number of records). Several studies have already demonstrated the crucial role of propagule pressure, and especially the number of new immigrants, on the colonization success of NNS (Cassey et al., 2018;Alzate et al., 2020). The higher abundance of invasive NNS in the belts could thus result in a higher propagule pressure inside PAs, and a subsequent higher probability of establishment of invasive NNS in PAs. Four NNS -R. rugosa, R. japonica, N. caerulescens and B. vulgaris -were all among the top 9 NNS in 5-km belts but not among the top 9 NNS in PAs. These are clear candidates for future colonization of PAs. This information is of relevance for managers of PAs to remain vigilant to future non-native colonizers.
In conclusion, our study strongly emphasizes the role of colonizing events from the surroundings of PAs in shaping NNS communities inside PAs. Both the abundance and the composition of the NNS communities around PAs influence NNS within PAs. Moreover, our study also reveals differences which are highly relevant for the conservation of PAs, such as the overrepresentation of invasive NNS within PAs. For all these reasons, we strongly suggest expanding the focus of NNS management within PAs to beyond PA boundaries as recommended by Monaco and Genovesi (2014). Considering the significance of the impact of invasive NNS in PAs (Hulme et al., 2014), efforts in monitoring and controlling invasive NNS are required from the PA management authorities, but also surrounding landowners. Similar advice has already been provided for PAs surrounded by high human population densities (Spear et al., 2013;Liu et al., 2020) -our study generalizes and reinforces it. The focus on NNS in the vicinity of PAs is of relevance for future conservation strategies, especially to prevent incursions of NNS with severe ecological impacts.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. This data can be found at: https://artskart.artsdatabanken.no and https:// www.protectedplanet.net.

AUTHOR CONTRIBUTIONS
KH did the analysis. KH and AC wrote the first draft of the manuscript. All authors commented, approved the manuscript, and conceived the study.

FUNDING
KH has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 766417.

ACKNOWLEDGMENTS
We thank the editors of Frontiers in Ecology and Evolution for this Special Issue. Further, we would like to thank two reviewers for their supportive and constructive comments on the manuscript. We also thank the Norwegian Biodiversity Information Centre (NBIC) for providing their database on alien species. The opinions given herein belong solely to the authors and do not represent the views or policies of IUCN.
Supplementary Figure 1 | Spearman rank correlation index of non-native species (NNS) richness in protected areas (PAs) and the four continuous explanatory variables: Mean accessibility of PAs, surface area of PAs, NNS richness in 5-km belts and mean human population density in PAs.
Supplementary Figure 2 | Years of first record of NNS in PAs and within a distance of 0-5, 5-10, and 10-20 km from the PA (belt zones). The numbers above the boxplot indicate the numbers of PAs and associated belts in which the NNS were present. Small letters (a,b) indicate elements that are significantly differentiated from each other according to Dunn's Test with "fdr" adjustment following a significant Kruskal-Wallis Test (KW).