Butterfly Populations Damaged By Increased Invasive Plants Richness

The monitoring scheme

Switzerland has a surface area of about 41,000 km² and encompasses a large diversity of environmental gradients as about 50% of its area is mountainous (elevations ranging from 190 m to 4634 m a.s.l.). To obtain a representative sample of the butterfly and plant diversity along these extended gradients, we used the Swiss Biodiversity Monitoring program dataset (BDM). BDM monitors the biodiversity of Switzerland since 2001 and is composed of 393 sites of 1 × 1 km size that are regularly distributed over Switzerland (Fig. S1.1 in Appendix S1). Within each of these BDM sites, vascular plants and butterflies were surveyed between 2007 and 2011 using standardized methods (i.e. transects of 2.5 km along paths and roads within the 1 km² plots). Depending on elevation, the sites were sampled once or twice for plants, and four to seven times for butterflies; annually between April 21st and September 21st, with a time interval of at least 14 days (Altermatt et al., 2008). A total of 1916 (native and exotic) plant species and 187 butterfly species were identified (Fig. S1.2 in Appendix S1).

The black list of invasive exotic plant species was obtained from the Swiss National Centre of Plant Data and Information (www.infoflora.ch) and contains 56 highly invasive and dominant plant species exotic to Switzerland (sensu Richardson et al., 2000) which show both a high spread potential, and a demonstrated negative impact on native biodiversity, human health and/or economy. For the following analyses, we selected in this list the 31 species that were terrestrial and present in the BDM plant inventory (see list of the 31 invaders in Appendix S1; black listed invader richness ranging from 0 to 15 species per site).

Functional traits and phylogeny

We used a set of functional traits describing the butterflies’ range size (as a proxy for habitat specialization), larval diet specificity, number of generation per year and morphological traits related to dispersal. As a coarse measure of the butterflies’ European distribution range size, we used the CLIMBER dataset (which covers all Europe and has a resolution of 50 × 50 km; Schweiger et al., 2014). Diet breadth of the larval stages was based on previously published field-observed interactions between food plant species and Lepidoptera (data from Ebert, 1991; Altermatt, 2010). Diet breadth of the larval stages was classified as monophagous (only one food plant species), strongly oligophagous (> 1 food plant species, but only one food plant genus), oligophagous (> 1 genus but only one food plant family) or polyphagous (> 1 food plant family). We used information on the number of generations per year (after Ebert, 1991), describing each species’ voltinism. Finally, we measured the thorax width and forewing length on life-sized photographs of all species, to calculate wing–load ratios (i.e. ratio of thorax width to forewing length). Higher wing–load ratios are usually associated with higher dispersal abilities (Dudley & Srygley, 1994; Turlure et al., 2010; but see Sekar, 2012).

The butterfly phylogeny is based on a molecular phylogenetic analysis of the concatenated mitochondrial gene cytochrome c oxidase I (1532 bp) and the nuclear gene elongation factor 1 alpha (1725 bp) of 85% of European butterfly species including all but two Swiss butterfly species (a total of 425 sequences; M. Wiemers & O.Schweiger, unpublished). We used the maximum likelihood method based on a general time-reversible model, starting the heuristic search with trees obtained by applying the neighbour-joining method to the matrix of pairwise distances estimated using the maximum composite likelihood approach. A discrete gamma distribution was used to model evolutionary rate differences among sites. The phylogenetic tree with the highest log likelihood was chosen. This analysis was conducted with mega5 (Tamura et al., 2011).

Environmental variables

To explore the relationships between plant and butterfly richness along environmental gradients, we considered a large set of variables encompassing: topo-climate, habitat type and habitat diversity at different spatial resolutions, as well as the nearest distance to each habitat type (see Table S1.1 in Appendix S1, for a detailed list of variables). The topo-climatic variables included annual sum of growing degree-days, mean annual temperature, mean annual potential evapotranspiration, mean annual moisture index, annual precipitation sum, number of frost days during the growing season, slope, aspect and elevation range per 1 km² pixel.

Statistical analyses

As a first step, we reduced the number of environmental variables by means of an ordination approach that has the advantage of providing uncorrelated axes of variation between sites in the subsequent regression analyses. We applied a principal component analysis (PCA) and used the first six PCA axes as representative of typical environmental gradients. These axes accumulated to explain 76% of the inter-site differences (Fig. S1.3 in Appendix S1). The first PCA axis notably represented low temperatures, high proportion of forest cover and long distances to urban and agricultural areas (Axis 1). The second PCA axis represented high habitat diversity and high proportion of urban areas (Axis 2). The third axis was correlated to high proportion of wetlands (Axis 3), while the fourth axis was related to high precipitation (Axis 4), the fifth axis to the number of frost days (Axis 5) and the sixth axis to the amount of solar radiation (Axis 6; Fig. S1.3 in Appendix S1). It can be noted that the proportion of urban cover, a variable important for both plants and butterflies, was well represented by the PCA axes as the six PCA axes together predict 81% of the urban cover variance in our dataset.

We aimed at estimating whether butterfly richness was affected by invasive plant richness, in addition to environmental variables and native plant richness. But as both native and invasive plant richness are also influenced by environmental variables, we adopted an approach avoiding spurious correlations between butterfly and invasive plant richness due to common correlations to environmental variables (for instance, if both depend on temperature, they are likely to appear correlated). To do so, we implemented three analysis steps: first, we modelled native plant richness as a function of environmental variables (i.e. six PCA axes). Second, we modelled invasive plant richness with the environment and the residuals from the native plant richness model as predictors (i.e. the influence of the native plant richness that was not due to the environment). Third, we modelled butterfly richness as a function of the environmental variables, the residuals from the native plant richness model, the residuals from the invasive plant richness model (i.e. the influences of the native and invasive plant richness that were not due to the environment) and the interaction terms between the residuals from the invasive plant model and the environmental variables. For all models, we allowed environmental variables to have both linear and quadratic terms.

For each regression, we fitted a generalized linear model (GLM; using Poisson or quasi-Poisson distributions). Full initial models were optimized by means of a stepwise backward–forward variable selection based on AIC scores optimization (the model with the lowest AIC score was retained; note that we used quasi-AIC for the quasi-Poisson models as suggested by Bolker, 2014) with the aim to retain only the statistically relevant variables. All final model residuals were checked to comply with the normality assumption. The richness of invasive plant species was thus identified as being correlated to the overall richness of butterfly species when retained in the final set of variables, with an estimated coefficient being significantly different from 0. The explanatory power of the GLMs was evaluated by means of the regression R² (based on the proportion of deviance explained and corrected by the number of degrees of freedom), and variable importance was estimated using a 99-fold randomization procedure for each variable separately (Strobl et al., 2007). The principle behind these randomizations was to mimic the absence of the variable in the model while maintaining the original degrees of freedom. The importance is then estimated as the average difference in prediction accuracy (i.e. regression R²) with and without permuting the target variable. Note that, we did not find spatial autocorrelation in the residuals of any of the fitted models (see mantel correlograms in Fig. S2.1 in Appendix S2).

To test whether individual butterfly species were particularly (positively or negatively) affected by local richness of exotic species, we then built logistic regression models for each butterfly species separately. In these models, the presence and absence of the species was explained by means of both the preselected environmental variables (i.e. 6 PCA axes) and the residuals of the native and invasive plant richness models (i.e. the influences of the native and invasive plant richness that were not due to the environment). We again applied a stepwise backward–forward variable selection procedure based on AIC scores (starting from a full model). For each butterfly species, the invasive plant species richness was considered to have a significant influence if it was retained in the stepwise-optimized model and had an estimated coefficient significantly different from 0 (with a P-value corrected for multiple testing across each butterfly species; Benjamini & Yekutieli, 2001). Note that for statistical robustness, we analysed only those butterfly species that had at least 30 occurrences in our dataset (115 species in total, i.e. 62% of all butterfly species in the inventory).

It can be noted that all mentioned regressions rely on the assumption that at the landscape scale, invasive plant species may influence butterfly species occurrence and richness, but are not affected by butterflies in return. This assumption is driven notably by the repeated evidence of an effect of specific invasive plants on native butterfly species (e.g. Tallamy & Shropshire, 2009; Jahner et al., 2011; Davis & Cipollini, 2014), but none (to the best of our knowledge) about the landscape scale effects of native butterflies on invasive plants. The latter may be found though at the local scale.

Then, we classified the butterfly species according to their response to invasive plant richness into three groups: (1) those that were positively correlated (i.e. eventually profiting from invasive plant richness), (2) those that were negatively correlated (i.e. eventually suffering from invasive plants richness) and (3) those that did not reveal correlation with invasive plant richness. We then explored whether these species groups differed in their functional and phylogenetic characteristics. Specifically, we used ANOVA to test for functional differences in (1) species’ range size in Switzerland (as estimated by the number of sites in which they occur), (2) species’ range size in Europe, (3) species’ larval diet specificity, (4) species’ voltinism and (5) species’ wing–load ratio. The patterns of phylogenetic relatedness among species belonging to each of these groups were tested by means of the mean phylogenetic distance between all pairs of species (MPD; Webb et al., 2002). To quantify whether the species within the groups were more clumped (clustered) or spread (overdispersed) across the phylogeny than expected by random grouping, we applied a null model randomizing the species’ position on the tips of the phylogeny (999 repetitions), and we calculated the standardized effect size (SES), hereafter called MPD_SES. MPD_SES can vary between 0 (completely clustered species) and 1 (completely overdispersed species), where 0.05 and 0.95 significance thresholds were applied.

Finally, to visualize the areas where the butterflies profiting or suffering from exotic species are located, we mapped both the prediction (from the individual models generated in the second step) and observation of these species. To map the individual model predictions, we transformed the continuous habitat suitability predictions into binary prediction with a threshold maximizing the proportion of correctly predicted presences (i.e. model sensitivity) and the proportion of correctly predicted absences (i.e. model specificity).

All statistical analyses were performed in R version 3.1.2 (R Development Core Team, 2014), and the packages were as follows: ade4 (Dray & Dufour, 2007), mass (Venables & Ripley, 2002), bbmle (Bolker & R Development Core Team, 2014), ecodist (Goslee & Urban, 2007), presenceabsence (Freeman & Moisen, 2008).

Plant and butterfly richness along the gradients

The selected set of environmental variables successfully predicted plant and butterfly richness, with explained deviances of R² = 0.44 for native plant richness, R² = 0.65 for invasive plant richness and R² = 0.75 for butterfly richness. These three groups all responded negatively to the number of frost days during the growing season (PCA Axis 5, Fig. 1; see Table S2.1 in Appendix S2 for the full list of estimated parameters of each model), but they responded differently to all other variables (Fig. 1).

The two most important variables for explaining butterfly richness in Switzerland were environment-independent native plant richness (40% with a positive coefficient; hereafter called native plant richness) and the annual number of frost days (Axis 5, 26%; Fig. 2). The environment-independent number of invasive plants (hereafter called invasive plants richness) nonetheless improved the butterfly model significantly (P-value < 0.001) and contributed to 4.1% of the explained deviance. The richness of invasive plant species revealed an overall negative effect on butterfly richness (P-value < 0.001), and significant interactions with Axis 1 (temperature, forest cover, distance to urban and agricultural area), Axis 3 (wetland cover) and Axis 5 (annual number of frost days) as shown in Fig. 2 and Table S2.1 in Appendix S2.

Characterizing butterflies that suffer from invasive plants

Among the 115 most abundant butterfly species modelled, we detected 28 species (24%; 16 genera) whose occurrence probability significantly decreased when invasive plant richness increases (Fig. 3b; see the full list of species in Appendix S2), but not a single species was positively affected. These butterfly species negatively associated with invasive plant richness had significantly smaller wing–load ratios than butterflies that were not affected by invasive plant species richness (P-value < 0.001; Fig. 3a). All other tested traits showed no significant differences between butterfly categories: species range sizes in Switzerland (P-value = 0.60) and in Europe (P-value = 0.82), larval diet specificity (P-value = 0.06) and average number of recorded generations per year (P-value = 0.16). The phylogenetic distance between butterflies that are negatively related to invasive plant richness was significantly lower than expected by chance (MPD_SES = 0.016). Note that we could not consider the two burnet moth species (Zygaenidae family; both negatively related to invasive plant richness) in this phylogenetic estimate because they were not included in the phylogeny.

The butterfly species that are negatively related to invasive plant richness were particularly present (both in number and as proportion of the local butterfly richness) in the mountainous regions (i.e. north-western Switzerland in the Jura mountains and southern Switzerland in the Swiss Alps; Fig. 4).

Invasive plant richness negatively correlates with butterfly richness

Invasive plant and butterfly species richness did not respond in the same way to the environment, with invasive plants preferring warmer and more urbanized sites compared to butterflies (Fig. 1, Axis 1). In fact, most invasive plants tend to originate from warm regions (Pyšek et al., 2003; Walther et al., 2007; Hulme, 2009) and are usually more abundant in urban and agricultural areas due to the high levels of disturbance in these areas and the high, human-induced propagule pressure (Chytry et al., 2008; Cabra Rivas et al., 2015; Gallien et al., 2015). Additionally, invasive plant richness was a negative and significant predictor of the overall butterfly species richness (even after removing the effect of the environment and native plant richness on invasive species richness).

We also found that the negative influence of invasive plant richness was specifically important in cold sites close to forested areas, in sites of low wetland cover, and in sites where the annual number of frost days is particularly high (see interactions in Fig. 2). These results indicate that in unsuitable conditions for most of the butterfly species, invasive plants may additionally decrease the local resources for the butterflies, for instance by replacing beneficial native plant partners, by attracting predators or by being toxic for the butterflies (e.g. Tallamy & Shropshire, 2009; Davis & Cipollini, 2014). However, in sites with a high proportion of wetland cover, the invasive species richness had locally a seemingly positive effect on butterfly richness. These sites, however, are generally only suitable for a small set of highly specialized butterflies (see Axis 3 in Fig. 1), which are generally rare. Yet, we could not confirm these specific positive effects of invasive plants on butterfly richness at the butterfly species level, as our analyses were limited to the more common butterflies, while wetland habitat specialists are usually rare.

Furthermore, we found that the measure of invasive plant richness is a better indicator of butterfly richness than other variables that are sometimes used, such as the level of habitat diversity, or the proportion of urban and wetland cover (e.g. Concepción et al., 2015; see Fig. 2). The explanatory power of invasive plant richness can, nonetheless, represent either a direct link to the butterfly physiology, with exotic species providing toxic resources for larvae and/or adult butterflies, or a missing covariate in the model, such as the frequency of disturbance in the sampled sites. It is essentially impossible to completely rule out the possible effect of missing covariates in empirical studies. However, we have used many covariates that are usually good proxies for disturbances (e.g. proportion of urban and agricultural fields, as well as habitat diversity). Also, the explanatory power of the model is high (R² = 0.75), and we found no spatial autocorrelation in the residuals (see Fig. S2.1 in Appendix S2). Therefore, it seems unlikely that an important covariate with spatial structure is missing.

It should be noted that we only had information available on the richness and not on specific abundances of invasive species per site (such as relative abundance or cover proportion). Substituting richness with specific dominance information in our analysis would likely reinforce our observed patterns (as improving the data quality should increase the explanatory power). Such information would be particularly useful to identify priority targets for conservation measures (e.g. if one invader was highly dominant, it would greatly influence the presence and distribution of butterflies and should thus become a target for invasive management plans).

Least mobile butterflies are most affected by invasive plants

Not all butterfly species responded in the same way to invasive plants. Specifically, 24% of the analysed butterfly species were less likely to occur in sites with a high number of invasive plants. These vulnerable butterflies are significantly less mobile than the butterflies that are insensitive to the presence of invasive plant richness. This suggests that more mobile species are less affected than less mobile species because the former might more easily reach suitable host plants and escape sites affected by invasive plants (see also Warren et al., 2001). The higher sensitivity of dispersal-limited butterflies to invasive plants furthermore suggests that the invasive plants do not only use local resources that are not used by the native plant species, but replaces native plant species on which the butterflies depend.

Butterflies’ range size, specialization and voltinism do not seem to be significant drivers of their sensitivity to invasive plants. Nonetheless, as the sensitivity of butterfly species to invasive plant richness is a characteristic that is clustered in the butterfly phylogeny, additional functional traits (that show phylogenetic signals) may also be affected by invasive plants (such as male mate-searching behaviour or female oviposition behaviour; Pavoine et al., 2014). Even though it may be difficult to quantitatively analyse many important butterfly traits that can be related to their susceptibility to invasive plant species, the position of butterfly species in the phylogeny may be used as a first indication for their general sensitivity. The phylogenetic clustering of sensitive butterfly species has important further implications, as it indicates that some phylogenetic clades – such as the genera Polyommatus and Argynnis – are specifically sensitive to increasing numbers of invasive plants. As a consequence, increased introduction and invasion of exotic plant species can lead to a decrease in the phylogenetic diversity of local butterfly assemblages.

Although we did not test for specific pairwise interactions between individual invasive plant and butterfly species due to our study design, our results can be used as a first screening to identify particularly vulnerable butterfly species that would ideally be studied in a more specific pairwise plant–butterfly interaction context. Additionally, our screening approach can be used to target the butterflies most vulnerable to invasive plant richness. These species should be prioritized in conservation planning. However, we point to the fact that the already rare and endangered butterflies, which are likely also very vulnerable to invasive plants, could not be assessed in this study due to data scarcity.

The detected impact of invasive plants on both butterfly functional characteristics and phylogenetic diversity might lead to functional and phylogenetic homogenization of butterfly assemblages with increasing amounts of exotic plant species. Functional homogenization of butterfly communities has been already documented across multiple functional groups (Eskildsen et al., 2015) and was found as a consequence of land use change (Ekroos et al., 2010). Similar results were also suggested for bird species assemblages in Europe in response to climate change (Thuiller et al., 2014). If the most sensitive butterfly species will go extinct, the associated loss of specific functional traits (or phylogenetic history) cannot be regained, which calls for careful conservation planning (Winter et al., 2013).

Our results also highlight the potential cascading effect of plant invasions on multiple trophic levels. Indeed, invasive plant spread does not only lead to a decline in the phylogenetic diversity of plants (e.g. Winter et al., 2009), but also to that of butterfly communities. This suggests that other herbivorous insect groups, such as beetles or flies, are likely affected similarly by invasive plants, which calls for repeated analyses of diverse insect groups and potential effects on higher trophic levels (Bezemer et al., 2014; van Hengstum et al., 2014; Litt et al., 2014).

To protect the Swiss butterfly diversity, it is important to both prevent the introduction and spread of invasive plants and to protect the areas most at risk for dispersal-limited butterflies. Targeting areas for butterfly protection can be achieved by prioritizing sites (1) with the largest number of butterfly species identified as suffering from invasive plants (e.g. the Central Alps of Switzerland), and (2) with the largest functional and phylogenetic diversities (e.g. using reserve selection algorithms; Arponen et al., 2005; but see Winter et al., 2013).