Linear and non-linear effects of goldenrod invasions on native pollinator and plant populations

The increased introduction of non-native species to habitats is a characteristic of globalisation. The impact of invading species on communities may be either linearly or non-linearly related to the invaders’ abundance in a habitat. However, non-linear relationships with a threshold point at which the community can no longer tolerate the invasive species without loss of ecosystem functions remains poorly studied. We selected 31 wet meadow sites that encompassed the entire coverage spectrum of invasive goldenrods, and surveyed the abundance and diversity of pollinating insects (bees, butterflies and hover flies) and native plants. The species richness of native plants decreased linearly with goldenrod cover, whereas the abundance and species richness of bees and butterflies decreased non-linearly with increasing goldenrod cover. However, no statistically significant changes across goldenrod cover were noted for the abundance and species richness of hover flies. Because of the non-linear response, goldenrod had no visible impact on bees and butterflies until it reached cover in a habitat of about 50% and 30–40%, respectively. Moreover, changes driven by goldenrod in the plant and pollinator communities were related to species loss rather than species replacement. We demonstrated that the impact of goldenrod cover on a habitat is not instantaneous. Habit management aimed at preventing the invasion process and alleviating its impact should take into account that, for the non-linear relationships, the negative impact can appear rapidly after crossing the threshold point.


Introduction
The introduction of non-native species to habitats increases consistently as a result of globalisation (Amano et al. 2016). Established alien organisms-invasive species-cause environmental changes that threaten native biodiversity and the human economy (Pejchar and Mooney 2009). The environmental changes induced by invasive species concerns the composition of invaded communities (Moroń et al. 2009). The effects may move through food webs to influence other trophic levels and different habitats (Simao et al. 2010;Lenda et al. 2018). In particular, the loss of plant species, primary producers, may impact key ecological processes such as ecosystem productivity , decomposition (van der Putten et al. 2013), or ecosystem resilience to disturbances (Richardson et al. 2007). As a result, the diversity decline at the level of the primary producers may impact associated arthropod species (Haddad et al. 2009). However, whether plants and arthropod communities change at all points along the invasion pathway (i.e. a linear response to invasion), or if there is rather a threshold beyond which the communities change as the invasive plant becomes dominant (non-linear response to invasion), remains unclear.
It is frequently assumed that the impact of an invading species is proportional (linearly related) to its density or abundance in a habitat (Yokomizo et al. 2009;Elgersma and Ehrenfeld 2011;Panetta and Gooden 2017). However, one may hypothesise that the density-dependent impact of invader populations may also elicit a non-linear response from native species (Crooks 2005;Panetta and Gooden 2017). Thus, there might be a threshold of density of invasive alien species at which negative effects on native communities appear (Fig. 1). It is expected that, if the threshold of density is high (Fig. 1), the number of species effectively endangered during the course of invasion tends to be underestimated and hence the consequences of invasions on biodiversity might be undervalued.
This study explores the pattern of the densitydependent impact of invaders using pollinating arthropods. Pollinators are key components of ecosystems providing various ecosystem services important for agriculture and the human economy (Mouquet et al. 2012). The effect of invasive plants on a pollinator community is rather equivocal (Bjerknes et al. 2007). The invasion of alien plants might affect pollinator populations by decreasing the diversity of native plants, and subsequently the pollinator food-base (Moroń et al. 2009). At the same time, alien flowering crops might boost pollinator population sizes by increasing their resource availability (Tepedino et al. 2008). Thus, a non-linear impact of invaders on native pollinators may be caused by the buffering characteristics of the environment (McCary et al. 2016). For example, plant-pollinator food webs may exhibit high functional redundancy (Memmott et al. 2006). Thus, even though an invasive plant at lower densities may eliminate several native plant species, high functional redundancy of the food web reduces the probability that this will lead to a decline in pollinators. Moreover, a novel invader might not substantially reduce foraging options for arthropods (McCary et al. 2016). Less food-specialised pollinators may be able to find their preferred resource while ignoring the presence of the plant invader, resulting in zero or minimal changes to a pollinator community in response to plant invasions at low densities. Additionally, the non-linear impact of invasive species can be generated by the neutral or partially positive effects of invaders on natives (Hejda and Pyšek 2006;Hulme et al. 2012), for example if an invader constitutes a novel, collateral food-source at low abundance and density (Salisbury et al. 2015). Regardless of the direction of the invasive species impact on pollinators, a full understanding will require knowing whether this impact may be linear or nonlinear as the invader becomes more abundant.
As a model invasive alien plant we used North American goldenrods (Solidago canadensis and S. gigantea), among the most invasive species in Europe and Asia (Weber 2001). Flowers of this species are rich in pollen and nectar (Kabuce 2006), and are visited by a spectrum of pollinating insects in their native areas (Pavek 2012) and in invaded ones (Kabuce 2006). Earlier findings showed that, in meadow habitats invaded by non-native goldenrods, there is a significant decline of native plant diversity (Moroń et al. 2009;Pal et al. 2015). The goldenrod invasions also cause communities of native pollinators (Moroń et al. 2009;Fenesi et al. 2015), ants (Lenda et al. 2013;Kajzer-Bonk et al. 2016), beetles (de Groot et al. 2007; Baranová et al. 2014) or birds ) to disappear.
We selected sites ranging from 0 to 100% of invasive goldenrod cover, and surveyed the abundance and diversity of the most important pollinator communities of the temperate zone (bees, butterflies and hover flies). For better mechanistic understanding of the potential linear or non-linear pollinators' response to the invasive goldenrods, we also recorded the number of native plant species, which are the main resources for pollinators (Blaauw and Isaacs 2014). We expected that there would be a non-linear response of pollinator communities to increasing goldenrod cover. We identified threshold values at which the negative impact increases for different groups of pollinators. Moreover, we tested whether the response of plant and pollinator assemblages reflects the degree of goldenrod cover. Specifically, we address the following questions: (1) Is there evidence of nonlinearity in the numerical response of native plants and pollinators as the cover of goldenrod increases? (2) If so, what is the threshold goldenrod cover at which the native plants and pollinators start to show a decrease in abundance and species richness? And, (3) Are plant and pollinator community changes caused by species replacement (i.e., turnover) or elimination (i.e., nestedness) as goldenrod cover increases?

Study area
The study was carried out in the meadow landscape located in the valley of the Vistula river near the city of Kraków, southern Poland (Fig. 2). Using present and historical vegetation data, we mapped the meadows (Kornaś and Medwecka-Kornaś 1974;Dubiel 1995;Pepkowska 2002;Skórka et al. 2007;Moroń et al. 2009) where soil is periodically saturated with water. The meadows have similar geological history, climate and soil properties (Kornaś and Medwecka-Kornaś 1974). All the studied meadows originated from wet meadows dominated by Molinia caerulea, Galium boreale and Sanguisorba officinalis with a number of other plants suitable for pollinators (Kornaś and Medwecka-Kornaś 1974;Moroń et al. 2008). After the collapse of communism, the wet meadows were maintained mostly by an irregular management scheme, with intervals of several years between mowing, burning or grazing (Skórka et al. 2007).
The rapid change in management allowed the goldenrod invasion (S. canadensis and S. gigantea) to begin at all sites basically at the same time, 20 years Fig. 1 Conceptual model of density-dependent impact of invasive species on native communities. With the increasing density of an invasive species, the richness of native communities decreases. Decrease of native communities could be nonlinear (community 1 and 2) or linear (community 3) dependent on the invaders' density. Arrows at the top indicate the threshold points for community 1 and 2 before this study began (which was confirmed by historical vegetation data). Differences in invasive goldenrod cover appear to be the result of local, random factors occurring at sites, e.g. occasional mowing, grazing or burning (personal observations). Moreover, because goldenrods are herbaceous plants, we can assume that the impact of coverage is not dependent on the temporal component of their lifecycle development stages (time taken for development from seedling to mature plant), unlike, for example, the case for trees (Pawson et al. 2010). From all located meadows, we selected 31 patches, ensuring that invasive goldenrod cover ranged from 0% up to 100% (43.06 ± 35.32%; mean ± SD). To calculate goldenrod cover, the study sites were carefully inspected and all patches of invasive goldenrods were mapped with the help of GPS. The study sites were at least 1.07 ± 0.18 km apart and were surrounded by arable fields, forests and human settlements.
To control for the confounding effects of potential spatial gradients (sites' size and distance to arable fields, human settlements, meadows and woodlands) for the investigated impact of goldenrod cover on pollinators, the wet meadows were selected in a relatively homogenous landscape. However, we applied QGIS 2.18.14 (QGIS Development Team 2018) software to calculate the spatial gradients of the study sites and used Spearman's test to ensure that the characteristics did not correlate with goldenrod cover. The sizes of the sites (6.12 ± 6.48 ha) were not correlated with goldenrod cover (r S = 0.136, p = 0.274). The goldenrod cover of selected sites did not correlate with the sites' distance to the closest arable field (r S = 0.167, p = 0.369; 473 ± 516 m), woodland (r S = -0.277, p = 0.132; 259 ± 294 m), human settlements (r S = 0.175, p = 0.347; 239 ± 168 m) or meadows (r S = 0.082, p = 0.660; 125 ± 69 m).

Surveys
A 200 m transect was established in the middle of each site (Fig. 2;Pollard and Yates 1993;Westphal et al. 2008). Bees and hover flies were sweep-netted along the transect on each site during four surveys, first in June, second in July, third in August and fourth in September (overall 124 transect-walks). During the transect walks on each site, the collectors walked slowly, making 500 sweeps to standardise the sweeping effort. Sweeps encompassed all flowering plants at transects. Butterflies were counted on each transect, on four occasions from June to September of 2010, and were captured only for identification purposes. The duration of a single transect walk for butterfly assessment lasted 20 min. The order in which the transects were sampled was random. Each transect Fig. 2 Location of 31 study sites in the Krakow region, Southern Poland was visited during different parts of the day (9:00 a.m. to 6.00 p.m.) throughout the season to cover entire period of pollinator activity during a day. The native plants were surveyed at five permanent plots (each with an area of 3.14 m 2 ), distributed at each transect with a distance of 50 m between the plots. The number of native plant species was surveyed twice during the study: first in the pre-blooming period (beginning of May) and second at the beginning of the goldenrod flowering period (in the middle of July).

Statistical analysis
To study the relationship between goldenrod cover and its impact on meadow pollinating insects and the potential role of native plant richness in shaping the effect, we used generalised additive mixed models (GAMM) with Poisson error and a logarithmic link function. Applying additive models, we analysed the expected non-linear relationships between the response and the explanatory variable (Wood 2006). We employed separate models to check whether the per site species richness and abundance of bees, butterflies and hover flies, and the richness of plant species, depended on goldenrod cover (seven models overall). Because the effect of the goldenrod cover on pollinators was consistent across the pre-flowering and flowering periods (''Appendix 1''), and as the effect on native plants was not dependent on whether plants were forbs or graminoids (''Appendix 2''), we pooled the data for each group (bees, butterflies, hover flies and plants). We fitted the goldenrod cover with the help of splines, to allow for the possible non-linear relationship between the variables and the cover. The possible spatial dependence of the data was addressed by including the interaction of longitude and latitude in GAMMs fitted with thin plate regression splines (Wood 2006). Using this procedure, part of the variation of a response variable is explained by the spatial location of a given transect, which makes the residuals from GAMMs spatially independent (Wood 2006). We based the parameter estimation on the ''mgcv'' (Wood 2006) computing models implemented in R (R Development Core Team 2016).
Moreover, to identify a threshold value of goldenrod cover where the slope of function described its impact on pollinators and plant changes we used estimation of regression models with break-points. The break-point values were calculated for all significant non-linear relationships revealed by GAMMs. To this point we applied the R (R Development Core Team 2016) package ''segmented'' (Muggeo 2015).
Differences between species assemblages can be an effect of species replacement between sites (turnover, which leads to different assemblages) and species loss from site to site (nestedness, which leads to the poorest assemblage being a strict subset of the richest one, Baselga and Orme 2012). To understand the mechanism of potential effects of goldenrod cover on pollinator and plant assemblages, we used the Jaccard dissimilarity index. This index takes into account the species that are different between a pair of sites, and is defined as the proportion of the overall pool of species present at each site. The index ranges from 0 for identical pool of species at all sites to 1 for no common species between sites. The dissimilarity indexes were computed for each pair of sites for pollinators and plants, and the total dissimilarity partitioned into two separate components accounting for the dissimilarity derived solely from turnover and that derived from nestedness (Baselga 2012). Then, with the help of partial redundancy analysis, we tested the goldenrod cover contribution to the dissimilarity derived from turnover and nestedness of pollinators and plants. The statistical analysis was performed using the packages ''betapart'' (Baselga et al. 2013)
Linear and non-linear responses of plant and pollinator communities were found. The number of native plant species decreases linearly with increasing goldenrod cover ( Fig. 3a; Table 2), whereas the he abundance and species richness of hover flies did not change across goldenrod cover (Fig. 3b, c, Table 2). However, the abundance of bees and their species richness decreased non-linearly with increasing goldenrod cover (Fig. 3d, e; Table 2). As with bees, a nonlinear decrease in abundance and species richness of butterflies occurred towards high goldenrod cover ( Fig. 3f, g, Table 2).
Threshold points were identified for non-linear responses of pollinators to goldenrod invasion. There was a negative effect of invasive goldenrods on bee species, and the number of individuals appeared at cover (± SE) of 50 ± 2% and 49 ± 1%, respectively (Fig. 3). A comparable detrimental impact of goldenrods on butterfly abundance appeared at cover of 31 ± 14% (Fig. 3). For the richness of butterfly species, the relationship, albeit non-linear, was more monotonic, with a threshold point at 43 ± 12%, (Fig. 3).

Discussion
An underlying assumption of the impact of invasive species on native ones is the proportional negative relationship between invader density and native population sizes (Dupont et al. 2009;Elgersma and Ehrenfeld 2011;Panetta and Gooden 2017). Such a linear relationship is more expected for direct interactions between invaders and natives. For example, low abundance of invasive plants can affect a native plant community via competition for space, soil resources and light (Crooks 2005). Accordingly, we found a negative, linear relationship between goldenrod cover and native plant species richness.
However, the expectation for indirect relationships, such as recorded in pollinating insects, is not so obvious . Indeed, we showed that the negative impact of invasive goldenrod cover on bee and butterfly communities is non-linear. This nonlinear response creates an threshold point and indicates that the invader effects accelerate only after a certain density threshold value is reached. We propose at least two food-web based mechanisms behind the effect of invasive plants resulting in the observed nonlinear impact. The first possible mechanism arises from the buffering potential of the environment such as food-web redundancy (Gilbert and Levine 2013). The second mechanism might be a result of the neutral or positive impact of the invaders on pollinators by increasing food-base (Hejda et al. 2009;Stout and Morales 2009). Both mechanisms may lead to a delayed abrupt collapse of populations (Kajzer-Bonk et al. 2016) or, alternatively, even reverse the negative impact of invasive plants on native pollinators. For example, patches with moderate cover of goldenrods can still provide sufficient key resources to maintain the lifecycle of native insects (Stout and Morales 2009). Pollinating insects, depending on their specialisation, may also switch between different plant species to find food and nesting sites. If flowering, the invasive plant can itself be an important food base for many pollinators (Bjerknes et al. 2007). We found a strong negative linear impact of goldenrod cover on the richness of native plant species, and a non-linear impact on pollinators. Thus, the buffering effect of the food-web might be a result of the food base redundancy of pollinators. Moreover, invasive goldenrods can be used as a food source by only some bee and butterfly species (Skórka et al. 2007;. Butterflies appear to be least resistant to the invasion, there was a negative effect on the number of individuals at sites with 30-40% goldenrod cover. Butterflies are herbivorous insects depending on different plant species during consecutive life stages. As adults, they feed on nectar of plants that may be different than plant species used by their often monophagous larvae (Brock 2015). Many butterfly species are recognised as food specialists for which food source redundancy is relatively low-70% of collected species were food specialised. Thus, bearing  ) and butterfly (f, g) abundance and species richness as predicted by GAMM models presented in Table 1. ns p [ 0.05, **p \ 0.01, ***p \ 0.001, 95% CI are marked with grey polygons. Red lines indicate estimated threshold points for non-linear relationships in mind the negative impact of goldenrods on native plants, and the dependency of butterfly life stages on different plant species, the impact of the invasive goldenrods on butterflies may multiply, leading to coextinction of host and butterfly (Koh 2004). Moreover, whereas some butterfly species adults can use invasive goldenrods as a nectar source, it seems that goldenrods are not used as food sources by their larvae (Ebert 2005;Buszko and Masłowski 2015). As such, the buffering mechanism seems to work only partially for butterflies, causing their relatively fast response to the invasion. Among other pollinators, bees have the closest, very often mutual, relationships with native flowering plants, as both larvae and adult bees are dependent on nectar and pollen collected from flowers (Brock 2015). Compared to butterflies, the number of bee species and individuals appeared to show a strong non-linear impact to goldenrod invasion. This relative robustness might indicate the strong food base redundancy of bees, as only 35% of collected bee species were food specialised. The negative impact at the level of 50% goldenrod cover may also indicate that bees utilise goldenrods, but as a supplementary resource rather than one complementary to the non-invaded meadows (Moroń et al. 2019). Moreover, for some foodspecialised bees, a patch dominated by one flowering species might be not preferred (Moroń et al. 2009).

Bees
Hover flies, among all the studied pollinator groups, are the least flower specialised because both adults and larvae feed on nectar and pollen as well as on other food sources (Brock 2015). However, the larvae of some species can be specialised, for example by requiring decaying wood to complete their lifecycle, or in the case of predators feeding on aphids (Brock 2015). Unlike bees and butterflies, which often have mutual relationships with a particular group of plants (Brock 2015), hover flies seems to be capable of using flowers of many species as food sources. The ability to use available resources might give an advantage, especially in habitats homogenised by few food sources. This may explain why, for hover flies, we found no direct link between the number of species and individuals, and goldenrod cover.
The pollinator and plant community assemblages seem to be affected as goldenrod cover increases. Changes in the community assemblages are a result of species loss rather than of the constitution of new communities by species replacement. Thus, in areas highly impacted by invasive goldenrods, only a part of the original pool of species can persist. The next step should be to identify traits which make some species less vulnerable to biological invasions.
Global environmental changes, along with the pressure of invasive species, has prompted a need to understand the trajectories of complex ecological systems (Mouquet et al. 2015) for better projections of biodiversity loss or changes and ecosystem functioning on the global scale (Barnosky et al. 2012). Unfortunately, biological invasions seem to be idiosyncratic, including non-linear relationships between invaders and natives, which makes prediction much harder (Richardson et al. 2000). However, if non-linear relationships are results of the understood intrinsic processes, the relationships can be incorporated into ecological predictions and plant management (Crooks 2005). The existence of the non-linear relationships also indicate some resistance of native communities to invasion, at least until the invasive species reaches a certain density. This prompts the recommendation that goldenrod cover below 30% for butterflies and about 50% for bees should sustain pollinator populations and their ecosystem services in meadow habitats under pressure of alien goldenrods. It also indicates that habitats and their pollinator communities, even when invaded, can be successfully restored to the point of relatively moderate goldenrod density, lowering costs of invader removal and the reintroduction of native species (Szymura et al. 2016).