Effects of competition and habitat heterogeneity on native‐exotic plant richness relationships across spatial scales

The biotic resistance hypothesis posits that greater native species richness limits invasions of exotic species. However, negative native‐exotic richness relationships (NERRs) may reverse with increasing spatial scale, seemingly refuting the hypothesis. Here, we explore the effects of species competitive interactions, environmental factors, habitat heterogeneity and vertical vegetation tiers on the NERRs across spatial scales in native forests.


| INTRODUC TI ON
Biological invasions are major threats to biodiversity with strong impacts on ecosystem functioning and the provision of ecosystem services (Vilà et al., 2011;Vilà & Hulme, 2017). Biotic resistance of ecosystems is widely considered a key factor influencing susceptibility to biological invasions (Iannone et al., 2016;Jeschke et al., 2012;Levine & D'Antonio, 1999;Nunez-Mir et al., 2017). This hypothesis, generally referred to as the "Biotic resistance hypothesis" or the "Biodiversity-invasibility hypothesis" posits that species-rich ecosystems are more resistant to biotic invasions than ecosystems with low biodiversity (Elton, 1958). Competitive interactions are thought to be among the main mechanisms involved in biotic resistance (Levine et al., 2004). The establishment of exotic plant species can be limited by native plant communities preventing access to essential resources such as light, nutrients and space. For instance, invasions of exotic plant species tend to be less common in undisturbed forests with a closed canopy and limited light availability in the understorey (Jagodziński et al., 2019;Jesson et al., 2000).
Other mechanisms such as direct and indirect negative effects of herbivory, predation and pathogens on exotic species establishment can also be associated with biotic resistance (Bufford et al., 2016;Eschtruth & Battles, 2009;Levine et al., 2004). The efficacy of biotic resistance mechanisms may vary depending on the characteristics of exotic species (Bufford et al., 2016;Nunez-Mir et al., 2017). For example, Martin et al. (2009) showed that although undisturbed forests were resistant to most invasions, they were susceptible to invasion by shade-tolerant plants. Greater native diversity can improve biotic resistance by maintaining the competitive pressure for diverse resources both spatially and temporally (i.e., niche complementarity) (Levine & D'Antonio, 1999) or by promoting herbivore diversity which increases the likelihood of an exotic plant species being eaten (Bufford et al., 2016;Levine et al., 2004).
Biotic resistance has been studied using a range of biodiversity measures, but the relationship between native and exotic species richness is the most commonly used (Jeschke et al., 2012). The strength and direction of the native-exotic richness relationship (NERR) were found to vary with the spatial scale (i.e., extent) of the studied area (Fridley et al., 2004;Shea & Chesson, 2002). At small spatial scales (e.g., within experimental plots of 10 × 10 m or less), NERRs are mainly negative, suggesting the presence of biotic resistance ( Figure 1). This appears to be driven by interspecific interactions such as competition for resources and space (Jagodziński et al., 2019;Tilman, 1997). However, some studies of local populations at small spatial scales have observed positive NERRs (Cleland et al., 2004;Stohlgren et al., 2003;Wiser & Allen, 2006). When the spatial scale of studies was extended, the NERR was most often positive (Davies et al., 2005;Herben et al., 2004;Iannone et al., 2016). This could be because external biotic and abiotic factors influencing both native and exotic species composition in parallel overwhelm the effects of species interactions (Davies et al., 2005;Fridley et al., 2007;Levine & D'Antonio, 1999;Naeem et al., 2000;Nunez-Mir et al., 2017). The reversal of the NERR from negative to positive with increasing scale was described as the "invasion paradox" by Fridley et al. (2007).
Many potentially confounding factors are involved in biotic invasions that can influence apparent biotic resistance outcomes across F I G U R E 1 Summary of hypotheses about factors explaining the invasion paradox. (a) The invasion paradox describes the reversal of the nativeexotic richness relationship from negative at small spatial scales to positive at large scales. (b) Biotic resistance may create a negative native-exotic richness relationship (NERR) at small spatial scale but has no effect at large scales. (c) Habitat heterogeneity can generate the positive NERR at large spatial scale while having no effect at small scales. The hypothesis is that (b) and (c) combine to create (a). The diagram is our summary of mechanisms given in Fridley et al. (2007) Exotic richness Small  were suggested to be explained by biotic acceptance and facilitation processes (Fridley et al., 2007;Nunez-Mir et al., 2017). Biotic acceptance may be observed when favourable habitats with greater resource availability allow the coexistence of high numbers of both native and exotic species (Stohlgren et al., 1999). Facilitation relationships between native and exotic species may offset competition in specific environmental conditions (Bruno et al., 2003). Factors that change with spatial scales such as habitat heterogeneity, propagule pressure and disturbances may also influence NERRs and lead to positive relationships with increasing spatial scale (Fridley et al., 2007) while not necessarily affecting local biotic resistance. Spatial heterogeneity may favour higher numbers of both native and exotic species by providing a wider range of ecological niches than found in more homogenous areas (Chesson, 2000;Davies et al., 2005). The level of habitat heterogeneity usually increases with the spatial scale of the observation, which by itself can generate a positive NERR ( Figure 1). High propagule pressure is known to increase exotic species establishments especially during early phases of invasions (Iannone et al., 2016;Wiser & Allen, 2006). Propagule pressure may explain a positive NERR at both small and large scales. However, at small scales, propagule pressure might have a weaker effect than biotic resistance allowing the observation of a negative NERR, whereas at larger spatial scales, more areas with strong propagule pressure may be included, overcoming biotic resistance and generating a positive NERR (Fridley et al., 2007). Natural or anthropogenic disturbances can reduce local diversity, reduce the competitiveness of native species and increase resource availability, all of which favour exotic species establishment (Eschtruth & Battles, 2009;Lembrechts et al., 2016). Increasing the spatial scale should increase the area of disturbed habitat which would result in higher invasion rates and, consequently, contribute to a positive NERR (Fridley et al., 2007). Two recent meta-analyses confirmed the scale dependency of the NERR which became more positive with increasing plot or grain size (Peng et al., 2019;Tomasetto et al., 2019). However, no negative NERR was detected in these meta-analyses, even at small spatial scales, probably due to the combined effects of studies that found a negative, neutral or positive NERR at small scales (the possible reasons for which we discuss below). Spatial scale was described as a combination of grain size and the spatial extent, and in both studies, grain size had a stronger effect on variation in the NERR than spatial extent. They also noted an unbalanced number of studies across regions with few located in the southern hemisphere.
New Zealand is among the most-invaded countries (Hulme, 2020), which provides an excellent opportunity to examine the NERR. More than 25,000 vascular plant species have been introduced, of which 2146 are naturalized with self-sustaining wild populations, almost equal to the number of native vascular plant species (2300) (Diez et al., 2009). Previous studies have investigated biotic resistance of native forest in New Zealand, mainly regarding single exotic species and at the local scale (Sullivan et al., 2005;Wiser & Allen, 2006;Wiser et al., 1998), leaving a gap for more general studies addressing multiple factors potentially involved in plant invasions. Several biotic and abiotic factors were found to favour exotic species establishment, such as anthropogenic or natural disturbance, low elevation, high soil fertility or propagule pressure (Jesson et al., 2000;Sullivan et al., 2005;Timmins & Williams, 1991;Wiser & Allen, 2006). However, intact native forest in New Zealand has been shown to be relatively resistant to invasion by exotic species (Jesson et al., 2000;Wiser & Allen, 2006). Interestingly, greater native diversity was found to favour exotic richness in native forest at a small scale (Wiser & Allen, 2006;Wiser et al., 1998). Both native and exotic species richness can be limited by light availability in old-growth forest with a closed canopy (Oberle et al., 2009;Timmins & Williams, 1991), but the creation of canopy gaps after disturbance (e.g., tree fall) may release both native and exotic species from competition, especially in the ground tier, and thus generate a positive NERR (Eschtruth & Battles, 2009;Iannone et al., 2014). Measuring forest biotic resistance by separating native tree species from other native species and controlling for forest type or localized disturbance may reveal a different pattern than the one observed by Wiser and Allen (2006).
Here, our aim was to test the biotic resistance hypothesis in New Zealand's native forests by assessing spatial scale effects on the NERR ( Figure 1) and to what extent this is affected by potentially confounding factors. We used the New Zealand National Vegetation Survey (NVS), Land Cover Database (LCDB) and climate data to examine the presence of the NERR and the invasion paradox in native forest plots and to test possible causes of the invasion paradox.
Specifically, we had three aims: • Aim 1: test whether focusing on native tree richness (rather than all native plant species) reveals competition effects driving part of the NERR and helps explain why some studies do not detect an NERR; • Aim 2: investigate the relationship between native tree and exotic plant richness across increasing spatial scales from 20 × 20 m plots to 128 × 128 km plot aggregations; • Aim 3: determine if considering potential confounding factors, such as 3a: propagule pressure, 3b: species competitive interactions and 3c: habitat heterogeneity, can explain a reversing NERR with increasing scale.

| Study area
Native forest covers about 23% of New Zealand's land area, ranging from the seacoast to the subalpine tree line (~1500 m) across diverse climatic conditions (Wardle, 1991). The composition of dominant tree species varies among regions with four main physiognomic groups or forest types being recognized: beech forest, beech-broadleaved forest, beech-broadleaved-podocarp forest and broadleavedpodocarp forest (Allen et al., 2013;Wiser et al., 2011). Undisturbed lowland forests are dominated by evergreen broadleaved tree species or conifers (mainly Podocarpaceae) (Wardle, 1991). Long-lived conifers (e.g., rimu, Dacrydium cupressinum; kahikatea, Dacrycarpus dacrydioides; kauri, Agathis australis) often occupy the upper canopy, which may be discontinuous, while broadleaved trees (evergreen angiosperms) may form a continuous canopy beneath that, with a sub-canopy and understorey of small trees, tree ferns, saplings and shrubs (Wardle, 1991). Evergreen southern beech forests (Nothofagaceae) dominate cooler and drier areas of New Zealand, mainly at higher elevations (Wardle, 1991;Wiser et al., 2011).

| Data sources and processing
We obtained native forest inventory data from New Zealand's National Vegetation Survey databank (NVS; www.landc are.cri.nz/ scien ce/nvs) (Wiser et al., 2001). In this study, we used a subset of For each plot, we calculated overall native and exotic species richness by combining all plant species present across all tiers.
"Native trees" (as opposed to shrubs, etc.) were defined as being present in Tier 3 or above, that is, ≥5 m in height, including vines and tree ferns but excluding epiphytes (Tier 7A) and used to calculate native tree richness for each plot. Exotic and native "ground tier" species richness was obtained from exotic or native species occurring at Tier 6 (0 to 0.3 m). As there were very few native forest plots with exotic species in tier 3 and above (16 plots out of 949 had 1 or 2 exotic tree/shrub species, and none had more than 2 tree species), we focused on exotic plants in the ground tier and their relationships with native plant species richness.
Native tree basal area and native species ground cover at Tier 6 were used as proxies for interspecific competition for light and space. Tree basal area was measured on stems of at least 2.5 cm diameter at breast height (Payton et al., 2004). Native ground cover was calculated for each plot, or group of plots in enlarged grid cells, by combining the cover classes using a modified Braun-Blanquet scale (1:<1%; 2: 1-5%; 3: 6-25%; 4: 26-50%; 5: 51-75% and 6: 76-100%) of each native species present at Tier 6. Plots described as native forest in the "Observed land cover" record (land cover observed during the survey), with at least one native tree species (≥5 m in height) and tree basal area information available, were kept for the analysis (N = 949).
To account for propagule pressure of exotic plants, we deter- To obtain meaningful groupings of LCDB5 classes, we created six new higher-level classes by combining multiple similar classes, resulting in the following: native woody vegetation, native grasslands, exotic woody vegetation, exotic grasslands, disturbed land and non-vegetated land (see Table S1). Exotic woody vegetation, exotic grasslands, disturbed land and non-vegetated land were consid-  (Wratt et al., 2006).
To study spatial scale effects on NERRs, we generated six spatial scales from the LUCAS grid. We created new enlarged grids by increasing the LUCAS grid dimensions from originally 8 × 8 km to 16 × 16 km, 32 × 32 km, 64 × 64 km and 128 × 128 km. Plots located within the enlarged grid cells were grouped together as a new spatial scale group (Table S2). We summarized the plot information by calculating the mean proportion of adjacent land cover, mean annual temperature, mean annual total rainfall, mean tree basal area and mean native ground cover for each new spatial scale group.
Species richness was calculated for all native plant, native tree, native ground species and exotic ground species grouping at each spatial scale using plant species present within an enlarged grid cell.
Heterogeneity of climatic conditions was used as a proxy for habitat heterogeneity. Total annual rainfall and mean annual temperature heterogeneity were calculated from the standard deviation of total annual rainfall and mean annual temperature across plots present within an enlarged grid cell.

| Statistical analysis
All models were run using Generalized Additive Models (GAM) with optimized smoothing parameters following the Restricted Maximum Likelihood (REML) method from the "mgcv" package (Wood, 2017) in R (R Development Core Team, 2018). GAMs were preferred over generalized linear models as they can account for linear and nonlinear relationships between variables. The suitability of a model was assessed by examining residuals, the number of basis dimensions for smoothed parameters using the "gam.check" function which performs diagnostics for a fitted GAM model, and the percentage of deviance explained (Wood, 2017). Collinearity between smoothed terms was checked for each GAM model by examining pairwise concurvity values. Pairwise concurvity provides a value between 0 (no collinearity) and 1 (high collinearity) between paired factors.
Explanatory variables showing a pairwise concurvity value above 0.8 were not included in the same model. GAM model selection was performed to obtain models with the best combination of variables relevant for each tested hypothesis using AIC and percentage of deviance explained (Table S3 and S4).
We first selected proxies for propagule pressure based on the most suitable combination of exotic or disturbed adjacent land cover types (i.e., exotic grassland, exotic woody vegetation, disturbed land and non-vegetated land) and radius distance (200, 1000 and 5000 m). A correlation matrix was generated using the "Hmisc" package (Harrell, 2021) to check for collinearity between the different predictors (Table S5). GAM models were first fitted with a Poisson distribution to predict exotic richness based on native richness, annual mean temperature, annual total rainfall and the proportion of exotic or disturbed adjacent land covers for each radius distance, but these showed signs of overdispersion ( Figure S2). Therefore, GAM models were fitted with a negative binomial distribution instead. The best combination of adjacent land cover types and radius distance showing the lowest AIC (Akaike Information Criterion) (Sakamoto et al., 1986) and homogenous range of values (e.g., disturbed land variables were either very small or large at small scales, increasing uncertainty about the sign of the relationship with exotic ground tier richness) was selected (Table S6).
Then, we examined the role of native tree richness in the NERR by comparing models predicting exotic ground tier richness from all native richness or native tree and native ground tier richness at the plot scale (Aim 1) and with increasing spatial scales (Aim 2). Annual mean temperature, annual total rainfall and the proportion of exotic adjacent land covers were included in the models as covariables to control for climatic conditions and propagule pressure. Native tree and native ground tier richness became highly correlated at larger spatial scales. Therefore, we decided to keep GAM models with these two variables separated at smaller scales only when investigating the effects of spatial scale on the NERR.
To analyse the potential confounding factors influencing the NERR at increasing spatial scales (Aim 3), we focused on the native tree-exotic richness relationship (NTERR). The influence of propagule pressure on the NTERR (Aim 3a) was investigated by including the proportion of exotic land cover to our models from plot to large spatial scale. The role of species competitive interactions (Aim 3b) was investigated across spatial scales by predicting exotic ground richness from native tree richness, tree basal area and native ground cover. The climate heterogeneity effect on NTERR (Aim 3c) was tested by including mean annual temperature and annual total rainfall heterogeneity as predictors at different spatial scales. Annual mean temperature and annual total rainfall were also included in all models to control for climatic conditions. To account for the difference in sampling effort, we included the number of plots per enlarged grid cell as a covariate.

| RE SULTS
Across the 949 plots, 1182 vascular plant species were recorded (including 120 identified only to genus) with 955 native and 207 exotic species (see Tables S7 and S8 for a list of most common species).

(c)
There was a marked difference in the NERR when native species across all tiers or only the tree tiers (>5 m) were used as the independent variable at the plot scale (Aim 1; Figure 2; Table S3 and S9). Although the NERR for all native plants (all plant forms in all tiers) versus exotics in the ground tier was statistically significant at p < .001 and positive (Figure 2a), the NERR based on native trees was negative (p < .001, Figure 2b). In addition, the NERR based on native plants in the ground tier only was positive (p < .001, Figure 2c), reflecting a positive correlation between native ground tier and exotic ground tier richness.
We observed a reversing slope of the NTERR (upper tiers above 5 m) with increasing scale (Aim 2; Figure 3a; Table S10). At the plot scale (i.e., when single 20 m 2 plots on an 8 km grid were assessed), there was a significant negative NTERR (p < .001). At an intermediate scale (32 km grid), the NTERR was statistically non-significant, and at the regional scale (128 km grid), it was positive (p = .03, Figure 3). When investigating the relationship between native and exotic ground tier richness, we observed a significant positive relationship across all spatial scales (p < .001, Figure S3, Table S11).
Hereafter, we focus on the NTERR in more detail to explore its reversal from negative to positive with increasing scale.
The amount of exotic grassland near plots had a significant positive effect on exotic ground tier richness across all spatial scales (Aim 3a; Figure 3b,c; Table S10). By contrast, an adjacent land cover of exotic woody vegetation had a positive effect on exotic ground tier richness only at a small scale (Figure 3c; Table S10). However, the pattern of the NTERR was still negative at the small spatial scale and positive at the large scale despite consideration of exotic adjacent land cover ( Figure 3). However, exotic ground tier richness showed a negative relationship with annual total rainfall at the small and intermediate scale (i.e., plot level and 32 km grid) and affected the NTERR. The positive NTERR at the large scale (128 km grid) became non-significant when annual total rainfall was included in the model (Table S10).
When controlling for competitive relationships between exotic and native species (tree basal area and native ground tier cover), the NTERR still reversed with increasing scale (Aim 3b; compare Figure 4a-c with Figure 3a, Tables S10 and S12).
However, the strength of the NTERR was weaker at the plot level when tree basal area and native ground tier cover were included in the model than without those terms (as in Figure 3) (Chi 2 reduced from 21.18 to 9.72; Tables S10 and S12 respectively).
Exotic ground tier richness showed a significant negative relationship with tree basal area, but only at small (p < .001) to intermediate scales (i.e., up to the 32 km grid, p < .001) (Figure 4d-f).
Similarly, a significant negative relationship was found between F I G U R E 3 Relationships at three different scales (8, 32 and 128 km grid) between log exotic plant richness in the ground tier (solid curve) and (a) native tree richness (for trees taller than 5 m), (b) the proportion of exotic grassland at 1 km radius, and (c) the proportion of exotic wood at 1 km radius. Hatched lines show 95% confidence intervals. Marks above the x-axis indicate the distribution of observations. Level of significance with non-significant (ns), *p < .05, **p < .001, ***p < .001 plot level, deviance explained went from 34% for the model with native tree richness, climatic conditions and exotic adjacent land cover types as predictors to 46.9% with native tree richness, climatic conditions, exotic adjacent cover types, tree basal area and native ground tier cover) (Tables S10 and S12).
Climate heterogeneity (i.e., mean annual temperature heterogeneity and annual total rainfall heterogeneity) influenced the NTERR across spatial scales (Aim 3c; Figure 5, Table S13). The negative NTERR at plot level became statistically non-significant at larger spatial scales if climate heterogeneity was included (Figure 5a-c showing NTERR from 16 to 128 km grid), whereas it was significantly positive at large spatial scales when climate heterogeneity was not considered (Figure 3c). However, when assessed individually, neither climate heterogeneity variable showed a significant relationship with exotic ground tier richness at large spatial scale (Figure 5f,i). Only mean annual temperature heterogeneity had to be kept in the model to observe a non-significant NTERR at a large spatial scale (128 km grid, Table S13). Interestingly, total rainfall heterogeneity showed a significant positive relationship with exotic ground tier richness at medium scale (i.e., 32 km grid, p < .001) which became statistically non-significant at the large scale (Figure 5d-f).

| Role of native tree richness
Our results show the importance of considering relationships between overstorey and understorey plant species when investigating biotic resistance in forests, as trees can influence exotic plant species invasion processes. Plots with a greater number of native tree species had lower exotic richness, despite the positive relationship between ground tier native and exotic richness as well as overall native and exotic richness. This latter positive NERR has been observed previously at the local scale in New Zealand native forest (Wiser & Allen, 2006;Wiser et al., 1998). In our study, our first key finding was the positive NERR being mainly explained by the positive relationship between native and exotic species found at ground tier. Native and exotic species in the understorey may respond to similar environmental conditions allowing both native and exotic species richness to be higher in more favourable areas (Shea & Chesson, 2002). As the models used in our analysis controlled for climatic conditions, tree richness and propagule pressure, other biotic or environmental factors could influence native and exotic richness such as tree species composition and forest type (e.g., native beech forest is known for its low native and exotic species richness) (Allen et al., 2013), soil characteristics (e.g., PH, fertility, moisture) (Gilbert & Lechowicz, 2005; Timmins & F I G U R E 4 Relationships at three different scales (8, 32 and 128 km grid) between log exotic plant richness in the ground tier (solid curve) and (a-c) native tree richness (for trees taller than 5 m) from models including tree basal area and native ground tier cover, (d-f) mean tree basal area (m 2 per ha), and (g-i) mean native ground cover. Hatched lines show ±95% confidence intervals. Marks above the x-axis indicate the distribution of observations. Non-significant (ns), *p < .01, **p < .001, ***p < .001 Our second key result was finding direct evidence that overstorey trees are the main contributor to forest biotic resistance through resource competition effects, as hypothesized by Iannone et al. (2016). Mature trees are known to have a direct impact on seedling and sapling survival near the ground, and hence on canopy replacement and exotic species invasions (Iannone et al., 2018;Martin et al., 2009). In older forests with mature trees, these fill the canopy space, limiting light availability in the understorey (Brockerhoff et al., 2003;Jagodziński et al., 2019;Oberle et al., 2009) and potentially monopolizing resources (e.g., water, nutrients) (Coomes & Allen, 2007), hence limiting exotic plant establishment and growth. The observed negative NTERR suggests that greater native tree diversity is a main driver of biotic resistance of native forest. High diversity of large native trees is likely to increase the occupancy of the canopy space and the exploitation of available resources (e.g., Sercu et al., 2017).
Furthermore, the exotic plants established in New Zealand are mainly associated with open and disturbed habitats as many were introduced as garden plants, for agriculture or as accidental introductions with seed imports (Hulme, 2020;Sullivan et al., 2005;Timmins & Williams, 1991). Although native tree diversity appears to enhance biotic resistance against most exotic species in New Zealand, some shade-tolerant exotic species can still invade forests (Martin et al., 2009) such as Mycelis muralis, which was the most commonly seen exotic species (Wiser & Allen, 2006).

| Invasion paradox
In this study, we confirmed the reversing NERR with increasing scale which has been described as the invasion paradox (Fridley et al., 2007). This relationship has previously been documented mainly considering richness of all native plants, mainly in nonforest habitats, and irrespective of relationships between overstorey and understorey tiers (Fridley et al., 2007;Herben et al., 2004;Tomasetto et al., 2019). This contrasts with our results which showed a reversing pattern driven almost entirely by the NTERR, as the NERR with all native plants was positive at all spatial scales.
Several different processes have been suggested as potentially influencing the NERR across spatial scales such as interspecific interactions (e.g., competition, facilitation), habitat heterogeneity, disturbance and propagule pressure (Fridley et al., 2007;Shea & Chesson, 2002). At a local scale, species competitive interactions apparently caused the negative NTERR, which is consistent with previous observations and with theory (Chen et al., 2010;Naeem et al., 2000). Both proxies for interspecific competition, tree basal area and native ground tier cover, had a negative relationship with exotic ground tier richness at small and intermediate scales. These findings underpin the importance of both the density of the tree F I G U R E 5 Relationships at three different scales (16, 32 and 128 km grid) between log exotic plant richness in the ground tier (solid curve) and (a-c) native tree richness (for trees taller than 5 m) in models also including rainfall and temperature heterogeneity, (d-f) rainfall heterogeneity (square root of total rainfall variance between plots from a similar group scale), and (g-i) temperature heterogeneity (km grid scale and square root of mean annual temperature variance). Hatched lines show 95% confidence intervals. Marks above the x-axis indicate the distribution of observations. Non-significant (ns), **p-value < .001 canopy and the density of ground cover in forest biotic resistance.
The combination of large native trees (Jagodziński et al., 2019;Martin et al., 2009) and dense native species at lower tiers may greatly reduce the availability of light and space, limiting the establishment of exotic species with typically low shade-tolerance (e.g., by leading to lower seed germination and growth) (Hulme, 2020;Jesson et al., 2000;Martin et al., 2009). Interestingly, exotic ground tier species richness showed a negative relationship with native ground tier cover, whereas the relationship with native ground tier richness was positive. Proxies for interspecific competition such as species cover could encompass mechanisms involved in biotic resistance other than species richness and provide different information about the level of resistance against exotic plant invasions of a studied habitat (Jeschke et al., 2012).
For instance, Iannone et al. (2018) investigated the contribution of trees to biotic resistance by focusing on maximum tree height, density and evolutionary relatedness with exotic cover as a proxy for exotic plant dominance, whereas exotic species richness was associated with species establishment.
At large scales, intraspecific competitive interactions were not influencing NTERR as the NTERR remained positive despite considering tree basal area. The detection of intraspecific interactions may be limited at larger scale by habitat heterogeneity and the size of the studied area, as more plant individuals not directly interacting are included (Fridley et al., 2007). Larger spatial scales can include a greater number of available resources, and heterogenous or disturbed habitats which may display low species competitive interactions (Davies et al., 2005;Fridley et al., 2007). Here, temperature heterogeneity and mean annual rainfall were influencing the NTERR at large spatial scales. This observation is consistent with previous studies suggesting that external factors and habitat heterogeneity, which affect both native and exotic species, could drive the positive NERR at larger spatial scales (Davies et al., 2005;Naeem et al., 2000;Shea & Chesson, 2002). The NTERR became clearly non-significant at the largest scale when temperature heterogeneity was considered, which supports theories that increasing habitat heterogeneity with increasing scale is a key element of the invasion paradox. This indicates that habitat heterogeneity may be a stronger driver of the positive NERR at large scale than mean environmental factors, which is consistent with the findings of Davies et al. (2005). Although Davies et al. (2005) used different environmental factors (i.e., soil characteristics, slope and aspect), they found that habitat heterogeneity strongly influenced the NERR at larger scale, not habitat favourability, which corresponds to the hypothesis that favourable habitat can support more native and exotic species. Previous New Zealand studies observed that fewer exotic species were found in high-rainfall areas (Ullmann et al., 1995). Wiser and Allen (2006) observed a negative latitudeexotic richness relationship in New Zealand native forest using a subset of the NVS. They suggested that the observed pattern was driven by the South Island where the west coast has higher rainfall and higher native forest cover, compared to the drier east coast which has high exotic grassland cover, a greater human population density and low native forest cover.
Propagule pressure, represented by adjacent exotic land cover types, had a strong effect on exotic richness at all spatial scales.
A greater proportion of adjacent exotic grassland cover favoured exotic ground tier richness at all scales from plot level to the largest scale (128 × 128 km). These results are consistent with the literature (Iannone et al., 2016;Timmins & Williams, 1991) and may support the hypothesis of Wiser and Allen (2006) which linked the negative longitude-exotic richness relationship with high exotic land cover on the east part of the South Island. However, we observed the invasion paradox pattern despite allowing for propagule pressure. Adjacent exotic vegetation provides propagules of exotic plant species explaining the positive relationship (Iannone et al., 2015;Wiser & Allen, 2006). However, pine forests adjacent to native forest in New Zealand were shown to have more native than exotic species in their understorey, especially in older stands (Brockerhoff et al., 2003;Forbes et al., 2019). Adjacent exotic woody vegetation may promote both exotic and native richness, but its effect may be overridden by exotic grassland with increasing scale as exotic grasslands cover greater areas (10.6 Mha) in New Zealand than exotic forests (2.1 Mha) (Brockerhoff et al., 2003;Stats NZ, 2022). These results reveal that propagule pressure is an essential factor when studying biotic resistance as it can be a confounding factor influencing exotic plant establishment (Iannone et al., 2015;Wiser & Allen, 2006). However, propagule pressure does not explain the invasion paradox as it does not vary with spatial scale. This indicates that, in our case, Fridley's et al. (2007) hypothesis that the NERR should be negative when controlling for extrinsic factors does not apply.
In conclusion, we were able to investigate processes involved in NERRs in New Zealand native forest that have not been previously teased apart. We tested several key hypotheses that have been proposed to explain the invasion paradox by using a single detailed and standardized dataset of vegetation plot data combined with environmental data which allowed us to build novel models including predictors not previously considered together. We found clear evidence that native tree diversity can favour forest resistance to exotic plant invasions in New Zealand at a local scale where interspecific competitive interactions take place. When considering temperature heterogeneity and mean annual rainfall, we did not observe the reversing NERR with increasing scale.
These results are consistent with the hypothesis that species competitive interactions drive the NERR at small spatial scales while habitat heterogeneity combined with key external factors such as rainfall drive this relationship at larger scales. Although the adjacent land cover type did not contribute significantly to the invasion paradox in our analyses, it is a relevant factor to include as a covariable. We also showed the importance of considering vegetation structure to test the biotic resistance hypothesis especially as it applies to forest ecosystems. Finally, our study supports previous findings that the invasion paradox can be considered as an accidental consequence of the study design, which inevitably has different factors influencing species richness at different spatial scales.