Secondary Invasions Hinder the Recovery of Native Communities after the Removal of Nonnative Pines Along a Precipitation Gradient in Patagonia

The removal of nonnative species can lead to re-invasion by nonnative species, especially in communities with multiple co-occurring invaders. Biotic and abiotic conditions shape community structure, reducing the predictability of nonnative management. We evaluated plant community recovery after the removal of nonnative pines with an emphasis on the effect of environmental conditions on the nonnative species response. We compared clearcuts (where pine plantations were removed), pine plantations, and native communities along a precipitation gradient in Patagonia. Nonnative richness and cover were higher in clearcuts compared to native communities along nearly the entire precipitation gradient, with the exception of the harshest sites. Compared to native communities, invasion resistance was lower in clearcuts in the wetter sites. Native richness and cover were lower in clearcuts relative to native communities along the gradient. Species composition in clearcuts diverged in similarity from native communities towards the wetter sites. Plantations showed an extremely lower richness and cover compared to both clearcuts and native communities. Our study highlights that clearcutting is an ineffective strategy to manage nonnatives aimed at restoring native communities and elucidates the importance of environmental context in management approaches. Taken together, our findings reinforce the important consideration of both the biotic and abiotic context of nonnative management.


Introduction
The management of nonnative species is a current challenge for ecological restoration [1][2][3], whose major goal is recovering the characteristics of an ecosystem that were prevalent before invasion, such as increasing biodiversity and restoring ecological functions [4]. The outcome of nonnative species management is highly unpredictable and the recovery of community structure and ecosystem functioning are hardly ever achieved or even evaluated [2,5]. An increasingly reported problem is that after the removal of a dominant nonnative species, other nonnative species invade the area, a process called secondary invasions [2,6,7]. Yet, most studies addressing the management of nonnative species focus on the management of a single-invader, without considering their community context [2].
The fact that many ecosystems are invaded by multiple co-occurring nonnative species, and that many factors can modulate nonnative species establishment, both contribute towards the unpredictability of nonnative species management [6]. Thus, a better understanding of the conditions that promote nonnative species invasion after the removal of a dominant nonnative can help predict management outcomes, as well as improve the allocation of management efforts [8].
Under certain conditions, the removal of nonnative species has led to successful outcomes [9,10]. However, it is not clear yet which conditions favored successful or unsuccessful restoration outcomes. Abiotic conditions can shape the variability in relative abundances and overall plant species composition in post-removal areas, which increases the unpredictability of nonnative species management [6]. Particularly, harsh environmental conditions (i.e., sites with extreme resource limitations where the occurrence of abiotic conditions that create rapid plant mortality is common, such as frost, extreme heat, and drought) may act as a strong filter for nonnative species. In fact, it has been found that harsh environments have a lower number of invaders than favorable environments [11][12][13][14][15]. For example, Sorte et al. [16] found that drought favored native over nonnative species. Secondary invasions in harsh sites will likely depend on nonnative species adapted to harsh conditions being present and able to respond rapidly under poor growing conditions [6,17]. Therefore, secondary invasions should take place less often under harsh environmental conditions than in more benign conditions. Biotic conditions may also shape the community response in different ways to environmental conditions. For instance, the diversity-invasibility hypothesis posits that more diverse communities exhibit greater resistance to invasions than less diverse communities [18,19]. More diverse communities have fewer unexploited resources reducing invasions via resource competition [20]. Further, resource competition may be stronger in more benign environments, consequently promoting biotic resistance [21]. Similarly, interactions among nonnative species can influence secondary invasions [7,22]. Competition among nonnative species can determine that the removal of a dominant nonnative releases other sub-dominant nonnatives from competition, thus favoring secondary invasions [7,22]. Additionally, an indirect positive interaction among nonnatives can drive the accumulation of nonnative species in the community, an interaction mediated by the reduction of native species abundance [22]. This positive interaction promotes secondary invasions after the removal of dominant nonnative species [6]. In this context, there is a need to develop general principles regarding invader interactions across varying environmental conditions so that secondary invasions can be anticipated and managers can allocate efforts toward pre-or post-removal actions [6].
Nonnative Pinaceae species (hereafter pines) have been planted in several regions of the southern hemisphere (e.g., New Zealand, Australia, South Africa, and South America) for forestry purposes and have subsequently invaded native habitats [23]. Both pine plantations and invasions produce a wide spectrum of changes in native ecosystems [23]. For example, pines have changed vegetation structure and fuel loads in Patagonian treeless ecosystems, which increase the intensity and frequency of fires [24]. As a consequence, changes in fire regimes reduce the recovery of nonnative species and promote further nonnative invasions [25,26]. Additionally, below-ground impacts may be more difficult to reverse, giving rise to both biotic and abiotic soil legacy effects (e.g., changes in soil nutrients, soil biota, or soil seed bank) that can drive changes in subsequent plant community structure and ecosystem processes [27,28].
Removal of nonnative pines (both planted and invasive) is a common management strategy around the world aimed at passively recovering native ecosystems, yet little is known about its efficacy [23]. For example, in Patagonia, many pine plantations are harvested for timber but not replanted owing to current bans on planting nonnative tree species. Moreover, many plantations are clearcut with the goal of restoration to native communities, especially in forest ecosystems [29]. However, it is well known that passive restoration to pre-existing states can be a challenge [30,31]. Removing nonnative pines can lead to undesired invasions of other nonnative species, halting the recovery of native ecosystems [30,32], and leaving vast areas with low timber productivity or conservation values. Therefore, assessing the effect of nonnative pine removal on plant community restoration, and understanding conditions that promote nonnative species in areas previously occupied by pines, is critical to properly managing pine plantations after timber harvest and ultimately restoring native communities.
The objective of our study was to evaluate plant community recovery after the removal of nonnative pine plantations and whether the effects of clearcutting varied with environmental conditions. We hypothesize that the previous presence of pines favors the establishment of nonnative over native species due to soil legacy effects (e.g., changes in soil nutrients, depletion of soil seed bank, changes in mycorrhizal communities) [23,30] and disturbance effects (e.g., an increase in resource availability, mainly light) [33]. We also hypothesize that the strength of these effects changes along a precipitation gradient, where they are weaker under harsher environments (i.e., drier areas) than in more benign environments (i.e., wetter areas) [6,14]. Additionally, we expect that steppe native species will better respond to clearcut conditions than forest native species, as light conditions in clearcuts are more similar to those of steppes than forests. Overall, we predict that secondary invasions should be higher and native community recovery lower in clearcuts in more benign sites. Since pine plantations are the prior state of clearcuts, we also evaluated nonnative species invasions and the similarity of plantations in comparison with clearcuts and native communities. These comparisons allow us to control for the effect of initial conditions (plantation understory) on clearcut community structure and to evaluate the impact of this land-use change on native communities, respectively.

Study Region
We conducted this study in Northwestern Patagonia, Argentina. This region is characterized by a steep west-east natural precipitation gradient caused by the rain shadow effect of the Andes, which acts as a barrier to the moist air coming from the Pacific Ocean [34]. Rainfall is concentrated between April and September and decreases from ca. 3000 to 500 mm per year over 100 km [35,36]. In this study, mean annual precipitation decreased from 1270 mm per year in the most mesic sites to 630 mm in the driest sites. Precipitation data for each site was obtained from Fetch Climate Web [37]. Mean annual temperature is 7.9 • C, with maximum temperatures occurring during January and February [35]. Vegetation shifts as mean annual precipitation decreases. Along this gradient, the wettest sites are temperate forests, dominated by Nothofagus spp. that are first replaced by Austrocedrus chilensis forest and matorral vegetation type in the forest-steppe ecotone, and finally by semi-arid grasslands or shrublands in the dry steppe ecosystem [38].
Three distinct physiognomic units occur from west to east: forests, shrublands, and steppes [34]. In the western area of the region, the Patagonian-Andean forest dominates; a vegetation unit dominated by deciduous, evergreen, and mixed forests [39]. The deciduous forest is mainly dominated by N. pumilio and N. Antarctica, which are restricted to the wettest and highest elevations of the gradient [39]. N. Antarctica also dominates stumpy forests in the driest and eastern part of the gradient [40]. Between 37.8 • and 47 • S, there are also forests of N. dombeyi, N. obliqua and A. chilensis [39], which are the most represented in the region spanned by our study sites. In the northern portion of this region, A. araucana appears as a subdominant species in these forests. The following trees and shrubs are also common: Lomatia hirsuta, Maytenus boaria, Schinus patagonicus, Azara microphylla, Aristotelia chilensis, Chusquea culeou, and Berberis sp. [39]. In the extra-Andean portion, shrubs increase and grasses decrease as mean annual precipitation decreases [34]. In this ecotone, we find the grained steppe that enters into the eastern sector of the deciduous forests, shaping a mosaic of both vegetation types. The vegetation cover is relatively high (64%) and it is dominated by Festuca pallescens and accompanied by Rytidosperma pictum, Lathyrus magellanicus, and some shrubs such as Senecio sericeonites and Azorella prolifera [39]. In the driest portion of the gradient, the typical vegetation is the grained-shrubby steppe where the typical vegetation is dominated by Pappostipa speciosa, Pappostipa humilis, Poa ligularis, and Poa lanuginosa, as well as by the shrubs Adesmia volckmannii and Berberis microphylla. This vegetation type has many variants according to the subdominant species [39].
In this region, many nonnatives have been introduced since European occupation. For example, several tree species from the Pinaceae family are spreading in the southern hemisphere, including Argentina [23]. In this region, Pinus contorta and Pseudotsuga menziesii are the main invaders of the native communities. Although conifers are naturally represented in these communities by two native trees (A. chilensis and A. araucana), there is no native species from the Pinaceae family in this region (e.g., all pines are nonnative). Besides pines, there is a high richness of nonnative species in the region [41,42]. Some species are only casual but others are highly invasive, such as Rosa rubiginosa and Cytisus scoparius [43,44].

Study Design
To evaluate if the effect of clearcutting on secondary plant invasion and community structure varied with environmental conditions, we surveyed 16 sites (Table S1) along a precipitation gradient in the 2016-2017 growing season. The mean distance between sites was~20 km. At each site, we selected three land-use types: (1) Clearcut: communities assembled after pine plantation removal. Clearcuts were considered to be different from others when previous pine species were different or when clearcuts had different ages. Clearcut age varied from two to eight years. It is well known that time since clearcutting is an important factor determining native vegetation recovery, and this could produce a bias in our results if there was a correlation between clearcut age and precipitation. We evaluated this and we did not find a correlation between the age of the clearcut and precipitation (r = −0.013, p-value = 0.96). Therefore, we did not find evidence of a possible bias in our results regarding clearcut age co-varying with precipitation; (2) Plantation: pine plantations that represented the situation previous to clearcut. All pine plantations surveyed were at a mature stage as our purpose was to represent the ecosystem state previous to clearcuts; (3) Native communities: areas dominated by native vegetation, with low levels of anthropogenic disturbance that represent a reference community. Within each land-use type, we randomly placed three observational plots (4 m 2 each) to assess plant community structure. In each plot, we recorded plant species composition and abundance (i.e., percent aerial cover per species) (Table S2). Species were classified by origin as native or nonnative following Zuloaga et al. [36].

Data Analyses
We evaluated the interactive effect of land-use type and precipitation on different descriptors of community structure: (1) native and nonnative species richness; (2) native and nonnative species cover; (3) proportion of nonnative species; (4) proportion of nonnative cover; and (5) Shannon diversity index based on species-specific foliar cover.
We tested the interactive effect of land-use type and precipitation by fitting a set of Bayesian hierarchical linear models. We modeled each community structure descriptor separately, and all models included land-use type and precipitation as predictors. To capture the hierarchical structure in the data (where plots were nested into sites), we set the land-use type variable as a categorical plot-level predictor and the precipitation variable as a continuous site-level predictor. While these models varied in their probability distributions, all of them were represented using similar deterministic functions that can be summarized as follows: Plot-level model: Site-level model: In Equation (1), α N represents the effect of native communities at each j precipitation level, whereas α C and α P are the analogous effect (i.e., effect size) of clearcut and plantation communities, respectively, compared to native communities (our reference community). As the model was fitted at two levels, parameters at the plot-level were allowed to vary with precipitation. Thus, β 1N , β 1C , and β 1P are the slopes of the site-level linear regression models (Equations (1.1)-(1.3)) for the native, clearcut, and plantation land-use type, respectively, and represent the rate at which the effect of land-use type changed with precipitation.
To describe variability around the above deterministic pattern, we used different probability distributions depending on the nature of the response variable (i.e., on the values that it can theoretically take). To model the species richness (count data), we assumed that the response variable drew a Poisson distribution [45,46]. The Poisson parameter (λ) was modeled as a linear function of community type and precipitation by means of a log link function (Code S1). To model species cover, we assumed that the response variable drew a Gaussian distribution, and modeled the parameter µ as a linear function of the same predictors (Code S2). To evaluate whether or not the richness and cover changed with species origin, we included this variable in the above models (Code S1 and S2). To model the nonnative richness and cover proportion (varying from 0 to 1), we assumed a Binomial distribution for the response variable [45,46]. We modeled the Binomial parameter ρ as a linear function of the predictors using a logit link function (Code S3 and S4). Finally, to model the Shannon diversity index, we assumed that the response variable drew a Gaussian distribution, where the parameter µ was a linear function of land-use and precipitation (Code S5). The response variables with continuous positive values (i.e., Shannon Index and cover) were modeled using Gaussian distributions as preliminary models employing log-normal distributions failed to converge. All these models were implemented in JAGS via the R package 'jagsUI'. We ran three chains with 10,000 iterations each discarding the first 5000 as burn-in.
To evaluate shifts in species' composition among land-use types, we performed Non-Metrical Multidimensional Scaling (NMDS). The ordination reduced the dimensionality of the distance matrix, and provided a first step for visualizing community dissimilarities [47]. We carried out a meta-NMDS from the R 'vegan' package [48] that generated an ordination of the Bray-Curtis distance matrix. Bray-Curtis distances represent how dissimilar two communities are, not only taking into account species composition (the list of species), but also the cover per species [47]. Bray Curtis distances were obtained with the 'vegan' package from a matrix in which the abundances per species registered at each of the three plots surveyed at each land-use type at each site were averaged. Additionally, we performed a permutational analysis of variance with the adonis function implemented in the 'vegan' package [48]. This non-parametric test allowed us to evaluate the interactive effect of precipitation and land-use type on the dissimilarity among communities. Finally, to evaluate if more invaded communities such as clearcuts tended to have less species turnover along the precipitation gradient, we estimated the Simpson beta-diversity index (B SIM ) for each land-use type along the gradient. If species invading clearcuts were the same along the precipitation gradient, we would expect that nonnative species composition would be less variable than total species composition along the gradient. For each land-use type, we estimated B SIM for total species composition and B SIM for nonnative species composition. Beta diversity indexes were estimated from the R package 'betapart' [49]. All analyses were conducted in R 3.4.3 (R Core Team, R Foundation for Statistical Computing, Vienna, Austria) [50].

Results
At the regional scale, we recorded 130 plant species (85 natives and 45 nonnatives) across all land-use types. We found 92 species (68 native and 24 nonnative species) in native communities, 83 species (46 native and 37 nonnative species) in clearcuts, and 31 species (19 native and 12 nonnative species) in pine plantations. The modeled native richness was lower in clearcut and plantation communities compared to the native communities at all precipitation levels ( Figure S1B,C). In contrast, our model showed that the effect of clearcutting and plantation on nonnative richness depended on the precipitation level (Figures 1D-F and S1E,F). Yet the modelling of the proportion of nonnative richness resulted in higher values in clearcuts and plantations than in native communities throughout the precipitation gradient (Figures 1G-I and S1H,I). In comparison with clearcuts, plantations harbored the lowest species richness, regardless of plant species' origin ( Figure 1B,C,E,F).  According to our models, species richness and proportion of nonnative species depended on the precipitation level in all land-use types, as reflected by the non-zero slopes in Figure 2. While native richness did not vary significantly with precipitation, nonnative richness increased as precipitation increased in native communities ( Figures 1A,D and 2A). This resulted in a lower proportion of nonnative richness in more benign (wetter) sites compared to harsher (drier) sites ( Figures 1G and  2B), as our model showed. In clearcut communities, both native and nonnative species increased towards more benign sites ( Figures 1B,E and 2A). Modeled nonnative richness was, on average, ~3fold higher in more benign sites (higher precipitation) than in harsher sites (lower precipitation) in this land-use type. Thus, the modeled nonnative richness was higher in clearcuts compared to native communities along nearly the entire precipitation gradient, with the exception of the drier sites ( Figures 1D,E and S1E). However, the proportion of nonnative species in clearcuts did not change with precipitation (slope close to zero) although the slope was marginally different compared to According to our models, species richness and proportion of nonnative species depended on the precipitation level in all land-use types, as reflected by the non-zero slopes in Figure 2. While native richness did not vary significantly with precipitation, nonnative richness increased as precipitation increased in native communities ( Figures 1A,D and 2A). This resulted in a lower proportion of nonnative richness in more benign (wetter) sites compared to harsher (drier) sites ( Figures 1G and 2B), as our model showed. In clearcut communities, both native and nonnative species increased towards more benign sites ( Figures 1B,E and 2A). Modeled nonnative richness was, on average,~3-fold higher in more benign sites (higher precipitation) than in harsher sites (lower precipitation) in this land-use type. Thus, the modeled nonnative richness was higher in clearcuts compared to native communities along nearly the entire precipitation gradient, with the exception of the drier sites ( Figures 1D,E and S1E). However, the proportion of nonnative species in clearcuts did not change with precipitation (slope close to zero) although the slope was marginally different compared to native communities ( Figures 1H and 2B). Our model indicated that, compared to native communities, plantations had fewer nonnative species in harsher sites but the difference in native and nonnative richness was diluted in more benign sites ( Figures 1D,F and S1F). plantations had fewer nonnative species in harsher sites but the difference in native and nonnative richness was diluted in more benign sites ( Figures 1D,F and S1F). Responses were considered different between native and nonnative species and among land-use types when 95 CI did not overlap with each other or with zero, respectively. Slopes in clearcut and plantation communities are relative to native communities (i.e., effect size). A positive slope indicates that species richness/proportion increased at a higher rate than in native communities, while a negative slope means that species richness/proportion decreased at a lower rate than native communities.
The modeled native plant cover was lower in clearcuts and plantations in comparison to native communities along the precipitation gradient ( Figures 3A-C and S2B,C). Instead, our model showed that nonnative cover increased with precipitation ( Figure A1A), following a similar pattern as nonnative richness. In the harsher sites, nonnative cover in clearcuts was similar to native communities, but in clearcuts, it tended to increase towards more benign sites ( Figures 3D,E and S2E). Proportion of nonnative cover decreased as sites became wetter in native and clearcut communities ( Figures 3G,H and A1B). In plantations, the modeled native and nonnative cover was close to zero along the precipitation gradient, but nonnative cover increased as sites became more benign compared to native communities ( Figures 3C,F and S2F). Responses were considered different between native and nonnative species and among land-use types when 95 CI did not overlap with each other or with zero, respectively. Slopes in clearcut and plantation communities are relative to native communities (i.e., effect size). A positive slope indicates that species richness/proportion increased at a higher rate than in native communities, while a negative slope means that species richness/proportion decreased at a lower rate than native communities.
The modeled native plant cover was lower in clearcuts and plantations in comparison to native communities along the precipitation gradient ( Figures 3A-C and S2B,C). Instead, our model showed that nonnative cover increased with precipitation ( Figure A1A), following a similar pattern as nonnative richness. In the harsher sites, nonnative cover in clearcuts was similar to native communities, but in clearcuts, it tended to increase towards more benign sites ( Figures 3D,E and S2E). Proportion of nonnative cover decreased as sites became wetter in native and clearcut communities ( Figures 3G,H and A1B). In plantations, the modeled native and nonnative cover was close to zero along the precipitation gradient, but nonnative cover increased as sites became more benign compared to native communities ( Figures 3C,F and S2F). Clearcut and native communities were significantly different in terms of species composition and species' relative abundance (F = 3.06, p-value = 0.001; Figure 4). The stress value obtained from the NMDS was 0.15, suggesting that the ordination was a good representation of the observed distances in the reduced dimensions. The differences among clearcuts and native communities were high (Bray Curtis distances higher than 0.5 in all cases) along the precipitation gradient. As the adonis test shows, the differences in community structure among land-use types varied with level of precipitation (F = 1.74, p-value = 0.001). Pairwise comparisons indicated greater differences between clearcuts and native communities in the more benign sites relative to the harsher sites (F = 1.83, pvalue = 0.006), as suggested by NMDS (Figure 4). Additionally, the diversity of clearcuts increased with precipitation (Figures 5 and A2); clearcuts were less diverse than native communities in the drier sites and became more diverse in the wetter sites ( Figures 5 and S3B). Likewise, the dissimilarity among plantations and native communities was also affected by the precipitation gradient (F = 1.80, p-value = 0.01). Finally, species turnover along the precipitation gradient (i.e., Simpson beta-diversity index (BSIM)) was high for overall species composition (BSIM-control = 0.876, BSIM-clearcut = 0.851, BSIM-plantation = 0.713) and for nonnative species composition (BSIM-control = 0.871, BSIM-clearcut = 0.814, BSIM-plantation = 0.692). In all land-use types, the identity of overall species composition and nonnative species alike varied along the gradient [40]. Clearcut and native communities were significantly different in terms of species composition and species' relative abundance (F = 3.06, p-value = 0.001; Figure 4). The stress value obtained from the NMDS was 0.15, suggesting that the ordination was a good representation of the observed distances in the reduced dimensions. The differences among clearcuts and native communities were high (Bray Curtis distances higher than 0.5 in all cases) along the precipitation gradient. As the adonis test shows, the differences in community structure among land-use types varied with level of precipitation (F = 1.74, p-value = 0.001). Pairwise comparisons indicated greater differences between clearcuts and native communities in the more benign sites relative to the harsher sites (F = 1.83, p-value = 0.006), as suggested by NMDS (Figure 4). Additionally, the diversity of clearcuts increased with precipitation ( Figures 5 and A2); clearcuts were less diverse than native communities in the drier sites and became more diverse in the wetter sites ( Figures 5 and S3B). Likewise, the dissimilarity among plantations and native communities was also affected by the precipitation gradient (F = 1.80, p-value = 0.01). Finally, species turnover along the precipitation gradient (i.e., Simpson beta-diversity index (B SIM )) was high for overall species composition (B SIM-control = 0.876, B SIM-clearcut = 0.851, B SIM-plantation = 0.713) and for nonnative species composition (B SIM-control = 0.871, B SIM-clearcut = 0.814, B SIM-plantation = 0.692). In all land-use types, the identity of overall species composition and nonnative species alike varied along the gradient [40].  Changes among land uses and along the gradient on the community descriptors were accompanied by changes in species dominance. Overall, clearcuts had greater relative cover of annual and perennial herbs and lower relative cover of shrubs and trees than native communities (Table A1). In clearcuts, the native tree Aristotelia chilensis and the non-natives Holcus lanatus (annual herb), Pseudotsuga menziesii (tree), Cirsium vulgare (annual/biannual herb), Pinus ponderosa, and Rumex acetosella had the highest cover values (Table S2). In contrast, in native communities, the most abundant species were the native perennial herb Chusquea culeou, the native trees Nothofagus antarctica, Schinus patagonicus, and Maytenus boaria, and the native shrubs Berberis microphylla and . Non-metric multidimensional scaling (NMDS) ordination plot of communities in two-dimensional scales. Each point represents the ordination score of a community, and the distance between any two points represents the difference between those two communities according to Bray Curtis distances. Communities that are closer together are more similar in composition, while communities that are farther apart are less similar. Ellipses represent 95% confidence intervals around the centroid of each land-use type. Colors indicate different land-use types: green for native communities, blue for clearcut communities, and red for plantation communities. Arrows point to the sites with higher precipitation. Point size indicates the mean annual precipitation of each site, where the larger the size of the point, the wetter the site.  Changes among land uses and along the gradient on the community descriptors were accompanied by changes in species dominance. Overall, clearcuts had greater relative cover of annual and perennial herbs and lower relative cover of shrubs and trees than native communities (Table A1). In clearcuts, the native tree Aristotelia chilensis and the non-natives Holcus lanatus (annual herb), Pseudotsuga menziesii (tree), Cirsium vulgare (annual/biannual herb), Pinus ponderosa, and Rumex acetosella had the highest cover values (Table S2). In contrast, in native communities, the most abundant species were the native perennial herb Chusquea culeou, the native trees Nothofagus antarctica, Schinus patagonicus, and Maytenus boaria, and the native shrubs Berberis microphylla and Changes among land uses and along the gradient on the community descriptors were accompanied by changes in species dominance. Overall, clearcuts had greater relative cover of annual and perennial herbs and lower relative cover of shrubs and trees than native communities (Table A1). In clearcuts, the native tree Aristotelia chilensis and the non-natives Holcus lanatus (annual herb), Pseudotsuga menziesii (tree), Cirsium vulgare (annual/biannual herb), Pinus ponderosa, and Rumex acetosella had the highest cover values (Table S2). In contrast, in native communities, the most abundant species were the native perennial herb Chusquea culeou, the native trees Nothofagus antarctica, Schinus patagonicus, and Maytenus boaria, and the native shrubs Berberis microphylla and Colletia hystrix (Table S2). Plantations showed lower total cover of both native and nonnative species but higher relative cover of trees than clearcuts and native communities (Table A1). The most abundant species in plantations were the non-native tree Pseudotsuga menziesii, and the rest of the species were notably less abundant than in the other community types (Table S2). Furthermore, we found a high level of turnover where few nonnative species occurred in more than half the sites (Table S2).

Discussion
Our results support the hypothesis that clearcut communities are more invaded by nonnative species than native communities. Proportion of nonnative richness and cover were higher in clearcut communities along the entire precipitation gradient. These results are similar to other cases previously reported where the removal of mature nonnative pines led to secondary invasions [10,30,32,51], but see [52]. However, Pauchard and Alaback [53] did not find high levels of invasions after pine clearcutting in the native range of the pine species. Disturbance generated by the removal of nonnative pines may increase resource availability and favor nonnative species, as suggested by the fluctuating resource hypotheses [54]. This occurs either by reducing resource uptake [33,54,55] or by increasing resource supply through residual biomass of the harvested trees [30,56]. However, these effects are more likely to occur immediately after clearcutting [56], which suggests that long-lasting legacy effects of pines may influence secondary invasions. In fact, pines produce below-ground impacts that can indirectly affect post-removal above-ground communities [57]. For instance, pines can reduce soil nutrient pools [58], decomposition rates [59], and soil pH [60,61]. Pines can also affect soil biota and native mutualisms [31,62,63]. Thus, ecological legacies of pines can indirectly promote the performance of nonnatives while hindering native species. For example, in New Zealand, the nonnative Pinus contorta altered biogeochemical cycles and increased ectomycorrhizal inoculum, which consequently generated a no-analog assemblage of species dominated by nonnative grasses and herbs after pine removal [30].
Native and nonnative species richness in clearcuts was higher at sites with the highest amount of precipitation. This suggests that native-rich communities tended to have more nonnative species than native-poor communities in clearcuts [12,13]. One possible explanation is that in drier sites, there would be fewer nonnative species adapted to the harsher conditions to be able to rapidly establish [6,14]. This would likely be due to introduction biases that altered nonnative species pools; nonnative species adapted to harsher conditions may be underrepresented in the nonnative species pool compared to nonnative species adapted to high-resource levels [64]. Instead, in the wettest sites, more benign conditions would not filter out the stress-tolerant species. Therefore, clearcut communities would have higher native and nonnative species richness and abundance [12]. The opposite occurred in the native communities where nonnative species richness decreased in wetter sites. This led us to hypothesize that under benign conditions, biotic resistance in native communities is higher than in clearcuts. In undisturbed conditions, an increase in biotic resistance in native communities may explain the lower nonnative cover and richness. Biotic resistance in native communities may be mediated by an increase in native cover. Previous studies have discussed the role of native species enhancing resistance to invasion owing to negative interactions among native and nonnative species [18][19][20]. Overall, young clearcuts showed greater nonnative invasions (had lower resistance to invasion) in more benign sites compared to native communities.
Regardless of secondary invasion, the success of passive restoration depends on the capacity of native species to (1) survive underneath pines and grow after removal or (2) recolonize the site from the soil seed bank or seed rain [9]. Here, we found that plantations had a negative effect on native richness and cover along the gradient, as has been previously reported in Patagonia [65,66]. This suggests that plantations are not harboring native species in either the harsh or in the more benign environments and that it is unlikely that they determine native community composition in clearcuts. However, under more benign environmental conditions than evaluated in our study, plantations can provide habitats for native species [67,68] and may accelerate passive restoration. For example, in more benign environments (1855 mm of mean annual precipitation) in New Zealand, Brockerhoff et al. [69] found a similar understory cover of native species in plantations and in the native forests. Seed banks allow new species that do not occur in the understory vegetation to occur after the removal of pines. Although our study did not directly address the role of seed banks influencing species composition, evidence suggests that soil seed banks are mainly dominated by nonnative species in pine-invaded communities around the world [10,70] and in other disturbed communities in Patagonia [71]. However, the importance of seed banks determining native vegetation dynamics in Patagonian communities is low [72,73] and variable along the precipitation gradient [73].
As we expected, young clearcut communities converge with native communities in the harshest sites. We hypothesized that the current high-light environment of the clearcuts favors the dominant native shade-intolerant species found in the harsher sites (steppe) and hinders native shade-tolerant species commonly found in the more benign sites (forest). Patterns of diversity suggest that the increase in native and nonnative richness in clearcut communities is driven by an increase in species abundance. Thus, the differences among clearcuts and native communities along the gradient could likely be explained by the higher proportion of nonnatives found in clearcuts across the gradient. Moreover, clearcuts presented a higher relative cover of annual herbs and lower relative cover of shrubs and trees in comparison with native communities, especially in the more benign sites. As native communities are mainly dominated by shrubs and trees, native species in clearcuts are notably different from those from native communities, especially in the wettest sites. We also found high species turnover along the gradient. Nonnative species in the harsher sites may not be a subgroup of those established in the more benign sites. Overall, one possibility is that clearcuts would need more time for passive restoration to succeed, particularly if the stage dominated by annual herbs is transient or facilitates the establishment of longer-live species typical of native communities. A more pessimistic scenario is that clearcut trajectories diverge from native communities, leading to alternative states, which can occur when plantations are burned [26]. Additionally, differences among plantation and native communities can be due to the extremely low diversity and understory cover found in plantations.
If the goal is to restore native communities following invasion, clearcut practices may not be an ideal technique to manage pines in Patagonia. Alternative practices, such as selective logging or techniques that leave dead pines standing (e.g., through girdling or herbicide application), may alleviate abiotic conditions (e.g., light, moisture) and promote native species compared to nonnative species if native propagules are not limiting [51,53,69]. Moreover, the management of current plantations (e.g., opening canopy) may increase understory biodiversity and accelerate clearcut restoration [65]. Additionally, management timing can influence restoration outcomes [23,66]. It has been found that the removal of pines allowed the regeneration of native communities in early stages of invasion. However, larger legacy effects appeared in later stages of invasion, hindering passive restoration [9,10,30]. In a more pessimistic scenario, return to the original native community may require additional interventions such as the re-introduction of locally extinct native species and their mutualists, or the modification of habitat conditions to make them more suitable for native species establishment [31]. It is important to note that secondary invasions may also generate economic problems by hindering the growth of desired species. For example, in its native range, Pseudotsuga menziesii (a species of interest for forestry worldwide) was negatively affected by the previous presence of the invasive nitrogen fixer Cytisus scoparius [74], a common nonnative species in clearcut communities in our study area.
Our work provides empirical evidence that furthers our understanding of the response of native and nonnative species composition to management under different environmental conditions. Thus, it may contribute to improve management approaches towards nonnative species and help the development of a theoretical framework for biological invasions [2]. Our results highlight the need to consider environmental context in the management of nonnative species. In harsh environments, a lower number of nonnative species respond positively to the removal of a primary invader compared to more benign environments, which also may have lower population and individual growth rates, as suggested by lower nonnative cover. Therefore, managers may have more time to manage or control the secondary invaders in harsher sites in comparison to more benign ones [2,6]. However, to achieve conclusive results, further studies should increase the sampling effort and time elapsed since nonnative removal, as well as conduct experimental studies that address the mechanisms underlying secondary invasions' patterns. Beyond that, based on our findings, it is not likely that clearcut communities that are already invaded by multiple nonnative species will recover and resemble native communities through passive restoration.

Conclusions
Removal of nonnative pines drives the secondary invasion of multiple nonnative species, altering plant species composition relative to those of native community assemblages. However, nonnative species richness and cover were higher in more benign (wetter) sites and clearcut communities were more similar to native communities in the harshest (drier) sites. The results of our work highlight the inefficacy of clearcutting to manage nonnative pines and restore native communities, especially in the wettest sites. Our conclusions also draw attention to the importance of the environmental context of management and reinforce recent arguments [2,6] that account for the biotic and abiotic context of nonnative species management as crucial.
Supplementary Materials: The following are available online at http://www.mdpi.com/1999-4907/9/7/394/s1, Table S1: Site locations (coordinates) with mean annual precipitation (mm) and time since clearcutting (yr) at each site, Table S2: Species cover (%) in 4 m 2 at each land-use type at each site; Figure S1: Effect size of each land-use type at each precipitation level for native richness, nonnative richness and proportion of nonnative species, Figure S2: Effect size of each land-use type at each precipitation level for native cover, nonnative cover and proportion of nonnative cover, Figure   Responses were considered different between native and nonnative species and among land-use types when 95 CI did not overlap with each other or with zero, respectively. Slopes in clearcut and plantation communities are relative to native communities (i.e., effect size). A positive slope means that cover/proportion increased at a higher rate than in native communities, while a negative slope means that cover/proportion decreased at a lower rate than native communities. Figure A2. Estimated rate of change of the Shannon diversity index along the precipitation gradient. Points represent the mean estimate of the slope of the hierarchical linear model that regressed Shannon with land-use type and precipitation. Vertical lines are credible intervals of 95% (95 CI) of the posterior distribution. Responses were considered to be different when 95 CI did not overlap with zero. For clearcut and plantation communities, the slopes represent the changes in the effect size of the treatment (i.e., relative to native communities) with precipitation. A positive slope means that Shannon diversity increased at a higher rate than in native communities, while a negative slope means that Shannon diversity decreased at a lower rate than native communities. Responses were considered different between native and nonnative species and among land-use types when 95 CI did not overlap with each other or with zero, respectively. Slopes in clearcut and plantation communities are relative to native communities (i.e., effect size). A positive slope means that cover/proportion increased at a higher rate than in native communities, while a negative slope means that cover/proportion decreased at a lower rate than native communities. Responses were considered different between native and nonnative species and among land-use types when 95 CI did not overlap with each other or with zero, respectively. Slopes in clearcut and plantation communities are relative to native communities (i.e., effect size). A positive slope means that cover/proportion increased at a higher rate than in native communities, while a negative slope means that cover/proportion decreased at a lower rate than native communities. Figure A2. Estimated rate of change of the Shannon diversity index along the precipitation gradient. Points represent the mean estimate of the slope of the hierarchical linear model that regressed Shannon with land-use type and precipitation. Vertical lines are credible intervals of 95% (95 CI) of the posterior distribution. Responses were considered to be different when 95 CI did not overlap with zero. For clearcut and plantation communities, the slopes represent the changes in the effect size of the treatment (i.e., relative to native communities) with precipitation. A positive slope means that Shannon diversity increased at a higher rate than in native communities, while a negative slope means that Shannon diversity decreased at a lower rate than native communities. Figure A2. Estimated rate of change of the Shannon diversity index along the precipitation gradient. Points represent the mean estimate of the slope of the hierarchical linear model that regressed Shannon with land-use type and precipitation. Vertical lines are credible intervals of 95% (95 CI) of the posterior distribution. Responses were considered to be different when 95 CI did not overlap with zero. For clearcut and plantation communities, the slopes represent the changes in the effect size of the treatment (i.e., relative to native communities) with precipitation. A positive slope means that Shannon diversity increased at a higher rate than in native communities, while a negative slope means that Shannon diversity decreased at a lower rate than native communities.