Global effects of non-native tree species on multiple ecosystem services.

Non-native tree (NNT) species have been transported worldwide to create or enhance services that are fundamental for human well-being, such as timber provision, erosion control or ornamental value; yet NNTs can also produce undesired effects, such as fire proneness or pollen allergenicity. Despite the variety of effects that NNTs have on multiple ecosystem services, a global quantitative assessment of their costs and benefits is still lacking. Such information is critical for decision-making, management and sustainable exploitation of NNTs. We present here a global assessment of NNT effects on the three main categories of ecosystem services, including regulating (RES), provisioning (PES) and cultural services (CES), and on an ecosystem disservice (EDS), i.e. pollen allergenicity. By searching the scientific literature, country forestry reports, and social media, we compiled a global data set of 1683 case studies from over 125 NNT species, covering 44 countries, all continents but Antarctica, and seven biomes. Using different meta-analysis techniques, we found that, while NNTs increase most RES (e.g. climate regulation, soil erosion control, fertility and formation), they decrease PES (e.g. NNTs contribute less than native trees to global timber provision). Also, they have different effects on CES (e.g. increase aesthetic values but decrease scientific interest), and no effect on the EDS considered. NNT effects on each ecosystem (dis)service showed a strong context dependency, varying across NNT types, biomes and socio-economic conditions. For instance, some RES are increased more by NNTs able to fix atmospheric nitrogen, and when the ecosystem is located in low-latitude biomes; some CES are increased more by NNTs in less-wealthy countries or in countries with higher gross domestic products. The effects of NNTs on several ecosystem (dis)services exhibited some synergies (e.g. among soil fertility, soil formation and climate regulation or between aesthetic values and pollen allergenicity), but also trade-offs (e.g. between fire regulation and soil erosion control). Our analyses provide a quantitative understanding of the complex synergies, trade-offs and context dependencies involved for the effects of NNTs that is essential for attaining a sustained provision of ecosystem services.

INTRODUCTION

Humans rely on the multiple services that ecosystems provide (MEA, 2005). Tree species play a key role in delivering ecosystem services, as they provide products (i.e. provisioning services, PES) such as timber, firewood, fibre, pulp or fodder. They also contribute to regulatory processes (regulating services, RES), such as climate regulation (via carbon uptake), soil formation and stabilization, and nutrient and water cycling (MacDicken, 2015; MacDicken et al., 2015). Trees also offer non‐material benefits (cultural services, CES), such as aesthetic or inspiration values; they are featured in the folklore, tales and legends of most human cultures, and contribute to people's sense of place (Kueffer & Kull, 2017; Mason et al., 2017).

With the increasing global demands of tree‐derived ecosystem services (MacDicken et al., 2015), many fast‐growing, stress‐tolerant or simply beautiful tree species have been extensively planted beyond their native ranges (Brundu & Richardson, 2016; Dickie et al., 2014a; Evans, 2009). Examples include non‐native conifers being planted worldwide for timber and pulp (Brundu & Richardson, 2016); Acacia, Eucalyptus or Pinus boosting land reclamation and sand dune stabilization worldwide (Evans, 2009; Griffin et al., 2011); legume trees (e.g. Acacia, Gleditsia or Prosopis) being used to prevent desertification and provide fodder or firewood in arid and impoverished regions of Africa and Asia (Shackleton et al., 2014; Witt, 2017); and shade and ornamental trees (e.g. Ailanthus, Jacaranda, Prosopis, Platanus and Robinia) providing amenities to rural and urban populations worldwide (Dickie et al., 2014a). Yet, many non‐native tree (NNT) species also contribute to landscape homogenization, reduce native biodiversity, and alter ecosystem functioning in undesirable ways (Cardinale et al., 2012; Gaertner, Richardson & Privett, 2011; Gamfeldt et al., 2013). NNTs may even result in direct negative impacts on human well‐being (i.e. ecosystem disservices, EDS) (Shackleton et al., 2016), such as toxicity or allergenicity. Therefore, NNTs also create threats to people's livelihoods and human well‐being, such as depletion of soil nutrients and water reserves (Castro‐Díez et al., 2012; Le Maitre et al., 1996; Shackleton et al., 2014), increased fire hazard (D'Antonio, 2000; Gaertner, Le Maitre & Esler, 2017; Nagler et al., 2005), damage to infrastructure and archaeological remains (Booy et al., 2017; Celesti‐Grapow & Blasi, 2004), or harm to human health (Nentwig, Mebs & Vilà, 2017; Schindler et al., 2015). These costs may be exacerbated when NNTs naturalize, and especially if they become invasive by spreading outside the areas where they were planted (Brundu & Richardson, 2016; Richardson & Rejmánek, 2011). Understanding how NNTs affect multiple (dis)services and how these effects correlate with each other is essential for enabling policy makers to balance costs and benefits, and manage conflicts over the use of NNTs (Dickie et al., 2014a; Pejchar & Mooney, 2009).

Most of our knowledge on the effects of NNTs on ecosystem services derives from local or regional studies, but these effects likely depend on the environmental conditions, history and cultural background of each region (Castro‐Díez et al., 2014a; Corbin & D'Antonio, 2011; Kueffer & Kull, 2017). Moreover, current knowledge is biased towards easy‐to‐study services and the most widely distributed NNTs (Hernandez‐Morcillo, Plieninger & Bieling, 2013; Hulme et al., 2013). Furthermore, many costs and benefits caused by NNTs are reported by different actors with multiple disciplinary backgrounds (e.g. foresters versus conservation biologists) using distinct publication types (e.g. scientific papers, reports or social media) (Krumm & Vítková, 2016). This paper presents the first global assessment of NNT effects on multiple ecosystem (dis)services, covering the three main categories of ecosystem services: regulating, provisioning and cultural services (de Bello et al., 2010; Haines‐Young & Potschin, 2013; MEA, 2005), and an ecosystem disservice. Specifically, we (i) evaluate the overall magnitude and direction of the effects that NNTs have on multiple ecosystem (dis)services, using native vegetation as a control; (ii) explore the role of environmental, geographical and socio‐economic factors as potential predictors of the variation of NNT effects on (dis)services; and (iii) identify synergies and trade‐offs among NNT effects on different ecosystem (dis)services. We address these questions through a global meta‐analysis of a comprehensive data set gathered from published scientific literature, country‐level economic and forestry reports, and social media.

MATERIALS AND METHODS

Given that different stakeholders may perceive the same service change as beneficial or detrimental, we avoid the value‐laden term ‘impact’ and use the more neutral one ‘effect’ to document objectively the changes produced by NNTs on ecosystem (dis)services (Jeschke et al., 2014). Thus, although throughout the manuscript we describe an effect as an increase/decrease of a (dis)service (Pyšek et al., 2012), we merely inform on the direction of the change, rather than judging the value of the change.

Data compilation

We focused on NNTs worldwide. We defined trees as ‘perennial woody plants with many secondary branches supported clear of the ground on a single main stem or trunk with clear apical dominance (including palms)’ (Richardson & Rejmánek, 2011, p. 789). ‘Non‐native’ trees were defined as tree species introduced (accidentally or intentionally) by humans to new geographic areas, considering the whole introduction–naturalization–invasion continuum (i.e. planted, naturalized and invasive species) (Richardson, Pyšek & Carlton, 2011). The ‘non‐native’ status of a tree was stated at the species level (i.e. ignoring the distribution of subspecies or varieties) and at the country scale, following regional and national floras and checklists (see online Supporting information, Appendix S1). For Brazil, USA and Canada however, non‐native status was considered at the state level due to their large size and environmental and biogeographic heterogeneity. Archaeophytes (i.e. species introduced before 1492) were excluded because in many cases they could not be clearly assigned to native tree (NT) or NNT at a country level. Hybrids between NTs and NNTs were also excluded for the same reason.

We covered the three most widely recognized categories of ecosystem services (de Bello et al., 2010; Haines‐Young & Potschin, 2013; MEA, 2005), i.e. regulating (RES), provisioning (PES) and cultural (CES) ecosystem services. We also considered an ecosystem disservice (EDS) – i.e. a negative impact on human well‐being (Shackleton et al., 2016) – specifically pollen allergenicity, which can be treated as a potential drawback from NNTs to human health (Vaz et al., 2017b). The selection of variables and data sources for each ecosystem (dis)service was achieved through an international participatory approach under the COST Action FP1403: Non‐native tree species for European forests – experiences, risks and opportunities (http://nnext.boku.ac.at), and relied on the possibility of worldwide coverage, cost and time efficiency, availability, and ease of communication across multiple countries. Data for RES were derived from an extensive scientific literature review, while data for PES, CES and EDS were collected from a thorough information search on reports and websites, among other sources (see Sections II.1a–d).

Regulating ecosystem services (RES)

We agreed on a list of 10 relevant RES that cover the benefits that people obtain from the capacity of ecosystems to regulate climate, floods, disease, waste, and water (de Bello et al., 2010; Haines‐Young & Potschin, 2013; MEA, 2005). For each RES, we identified a set of underlying target variables which covered ecosystem properties, processes, and/or traits of dominant species that underpin the capacity of ecosystems to regulate processes and mitigate effects of disturbances (de Bello et al., 2010; Quetier et al., 2007) (Table 1). We performed a literature search of scientific publications, using ISI Web of Knowledge (http://webofknowledge.com/) and Scopus (https://www.scopus.com). The search was carried out in December 2015 and updated in December 2016, covering the period between 1904 and 2015 with no restriction on language. Our search string included: (Exotic* OR Alien* OR invas*) AND (*native*) AND (tree*) AND the set of key words shown in Table 1 combined with ‘OR’.

Table 1

Target variables used as proxies for different regulating ecosystem services (RES). Variables include quantifications of ecosystem processes, ecosystem or community properties and traits of dominant plant species. The positive or negative sign beside each variable indicates the relation with the ecosystem service. The last column shows the list of specific key words used in the search in ISI Web of Knowledge and Scopus

	Target variables
Regulating ecosystem services	Ecosystem processes	Ecosystem/community properties	Plant species traits	Key words used in the literature search†
Climate regulation	Carbon sequestration + Biomass production+	Aboveground plant mass/C + Root mass + Soil carbon + Total plant mass + Tree basal area +	Chlorophyll concentration + Photosynthetic rate + Relative growth rate + Tree height + Trunk area/diameter + Trunk diameter increment +	Carbon sequestration, Carbon storage, Primary production, RGR, Growth rate, Photosynthetic rate, Chlorophyll concentration, Microclimate, Climate regulation, Canopy temperature, Wind
Fire‐risk prevention		Canopy fuel continuity – Canopy water content + Litter mass/depth – Litter water content + Understorey biomass –	Calorific value – Effective heat of combustion – Leaf moisture + Volatile compounds –	Fire, Fire frequency, Fire susceptibility, Fire intensity, Burning temperature, Fire spread, Forest fire, Wildfire Fire regime, Fire behavio?r, Fuel propert*, Flammability
Flood regulation		Flood frequency – Stream water velocity –		Flood frequency, Flood* Water velocity, River flow, Run?off, Flood protection, Flood defence, Flood storage, Flood generation, Flood detention, Flood event
Plague control		Forest plague frequency – Abundance insectivorous species +	Leaf lignin content + Polyphenol content +	Plague frequency, Disease frequency, Tree pathogens, Natural pest control, Pest control, Biological control, Biological pest control
Pollination			Pollinator visitation rate to flowers +	Pollinat*, Pollination service, Pollinator efficiency, Flower visitor, Zoophilous
Pollution regulation	Soil NOx emissions –		Concentration of heavy metals in tissues + Plant isoprene emissions – Plant monoterpene emissions – Plant NOx emissions –	Air purification, Air clean*, Pollut*, Contamination, Noise, BVOC, Biogenic emission*, Volatil*, Water quality, Water purification, Water clean*, Sequestration, Mining
Soil erosion control	Leaf litter production +	Litter layer mass/depth + Root mass per unit soil area + Understorey biomass +	Root depth +	Soil erosion, Weathering, Soil loss, Sediment, Root depth, Root density, Erosion protection, Soil stability, Sand stability, Root depth, Root density, Soil erodibility, Soil floor
Soil fertility		Canopy nutrient content + Carbon exchange capacity + Soil base saturation + Soil nutrient content +	Leaf nutrient content + Litter nutrient content +	N fixation, (Soil, Leaf, Leaves, Litter) AND (Nutrient*, Nitrogen, Phosphorus, CEC)
Soil formation	Infiltration rate of nutrients + Litter accumulation rate + Litter decomposition rate + Mineralization rate + Nutrient input by litterfall + Soil microbial activity + Soil respiration rate +	Litter layer mass/depth + Root mass per unit soil area + Soil invertebrate abundance + Soil organic matter +	Hyphal length + Litter C:N – Litter lignin – Litter lignin:N –	LMA, SLA, SLM, Nitrogen, Phosphorus, Lignin, Litter decomposition, Litter C:N, Litter C/N, Mineralization, Nitrification, Ammonification, Soil respiration, Microbial biomass, Soil organic matter, Soil compaction, RGR, Growth rate, Litter layer, Litter?fall, Soil invertebrates, Root specific length, Hyphal length
Water regulation	Canopy interception of rainfall + Evapotranspiration + Infiltration rate +	Canopy water content + Leaf area index + Litter layer mass/depth + Soil moisture + Water repellency –	Sap flow rate + Stomatal conductance +/– Transpiration rate + Tree water consumption rate + Water use efficiency +	Canopy water content, Soil moisture, Runoff, LAI, Litter layer, Evapotranspiration, Infiltration, Water recharge, Transpiration, Sap flow, Stomatal conductance, Water use efficiency

†Compound key words were introduced between inverted commas.

Publications retrieved from our searches were filtered according to the following criteria: (i) the study compares any target variable between a site dominated by a NNT and a control site with native vegetation; (ii) the NNT site and the control site are close to each other and have similar environmental conditions, according to the authors of the papers; and (iii) the paper provides means, deviation and sample size of the target variable. The final number of selected papers was 135 (Appendix S2).

Using data obtained from the selected publications, individual data sets were created for each RES (see Table 1). Each data set consisted of a table where rows were case studies and columns covered the mean, deviation and replication of a target variable in a NNT site and in a control site, as well as other explanatory variables, which were subsequently used to explain the variability of effect sizes (see Section II.1e). When the same paper provided measures on more than one target variable, we adopted any of the following decisions: (i) for variables associated with different ecosystem services, we kept them as independent case studies because each service was independently analysed; and (ii) for variables associated with the same service, we selected the one most directly related to the service or we aggregated all values into one to avoid pseudo‐replication (see Appendix S3 for details). In total, our data set contained 1688 case studies.

For each case study we calculated a standardized effect size using Hedges' d (Rosenberg, Adams & Gurevitch, 2000), i.e. the difference of mean values between the NNT and the control site, weighted by the number of replications and by the inverse of the variance (Appendix S4A). When the variable had a negative relation with the ecosystem service, we changed the sign of the effect size. A positive effect size thus indicates that the NNT site has a higher contribution to the particular ecosystem service than the control site and vice versa.

Provisioning ecosystem services (PES)

We considered two main categories of products obtained from areas where land cover was classified as ‘forest’: timber products (any kind of goods obtained from harvesting trees) and non‐timber products (any biological resource in woodland except timber) (MEA, 2005). Details on the types of products included in each category and specific sources of information for each country are shown in Appendix S5.

For timber products, we performed an intensive search of country reports where information on provisioning of harvested timber (m3 y–1) was available per species, at country or sub‐country (state or administrative region) level, for the period 2007–2015. The difference between the proportion of timber provided by NNTs and by NTs was calculated as the effect size for each country/region (i.e. case study). A positive value indicates a higher timber provision by NNTs than by NTs, and vice versa (see Appendix S5 for further details).

For non‐timber products, we considered 16 categories, following the classification of the Secretariat of the Convention on Biological Diversity (https://www.cbd.int/doc/publications/cbd-ts-06.pdf; Appendix S5). Fruits obtained from tree crops were excluded because they come from non‐forest land cover zones. Given that production of different non‐timber products is often species specific and may differ in orders of magnitude across products (e.g. cork production is expressed in t ha–1 y–1, while honey and edible fungi in kg ha–1 y–1) they could not be aggregated across species. Moreover, information on productivity per species at the region/state level was generally not available for several of these 16 categories. To overcome these limitations, we compared the uniqueness in the production of non‐timber products between NNTs and NTs. Uniqueness was recorded as the number of categories of non‐timber products (N = 0–16) that could be obtained only from NNTs or only from NTs at country level (e.g. in Spain and Portugal, cork is exclusively produced by a NT, so this product counts as unique to NTs; by contrast, honey can derive from both NTs and NNTs, thus not counting as unique to any group). We compiled data for a total of 16 countries, each one representing a case study.

To assess whether NNTs have higher non‐timber uniqueness than NTs in each country, we calculated the effect size using the log odds ratio under Peto's method (logOR ), which is frequently used in meta‐analyses to aggregate data reported as frequencies. The logOR is the difference between an observed value and the value expected by chance in a contingency table (Borenstein et al., 2009). In this case, the observed value was the number of non‐timber categories unique to NNTs, while the expected value was calculated from the number of NNT and NT species present in the country, expecting that both groups have the same uniqueness for non‐timber production (for further details on logOR calculation see Appendix S4B). A positive logOR value indicates that NNTs supply more unique non‐timber product categories than expected according to the proportion of NNT species in the country. Specific sources of information are shown in Appendix S5.

Cultural ecosystem services (CES)

CES cover non‐material benefits people obtain from ecosystems through spiritual enrichment, cognitive development, reflection, recreation, and aesthetic experiences (MEA, 2005). These benefits were grouped under five categories: recreation and ecotourism, aesthetic, inspiration, cultural heritage and scientific interest (MEA, 2005). For each CES category we selected at least one representative source of quantitative information (Vaz et al., 2018). For example, as sources for aesthetic value, we selected catalogues of ornamental plant dealers and catalogues of species present in urban parks (see Table 2 for all CES sources and Appendix S6 for further details). For each case study (e.g. each catalogue), we calculated the effect size using the logOR (as in Section II.1b). In this case, the observed value was the frequency of NNT species in the CES source; the expected value was calculated assuming that both NNT and NT species have the same chances of being included in the source (Appendix S4B). LogOR values higher or lower than 0, respectively, indicate an increase or decrease in the particular CES caused by NNTs. Data on CES were collected for different regions/states for most countries. For USA, only data for Hawaii could be collected. In total, our database contained 938 case studies from 13 countries.

Table 2

Sources of information associated with different categories of cultural ecosystem services (CES) and one disservice (EDS), pollen allergenicity. Calculation of the log odds ratio under Peto's method (logOR) was based on the difference between the observed non‐native tree (NNT) value in a given source (A) and the expected NNT value under the assumption that both NNTs and native trees (NTs) have the same chances of being included in the source: (A+B)×(A+C)/(A+B+C+D). For further details see Appendices S4B and S6. Information was collected mostly at the sub‐country level (state or administrative region). No., number.

			Values observed in the source associated with a given (dis)service		Values used as control
Cultural ecosystem (dis)services	Sources of information	Rationale	NNT in the source (A)	NT in the source (B)	NNT in the control (C)	NT in the control (D)
Aesthetics	Catalogues of ornamental plant dealers	Tree species offered by plant dealers are appreciated mostly for their aesthetic values	No. of NNT species offered in catalogues	No. of NT species offered in catalogues	No. of NNT species present in the country	No. of NT species present in the country
	Tree inventories of urban parks	Tree species exhibited in urban parks are included mostly for their aesthetic values	No. of NNT species present in inventories	No. of NT species present in inventories	No. of NNT species present in the country	No. of NT species present in the country
Recreation and ecotourism	Official tourism websites	Photographs from tourism websites were selected for the potential of NTs or NNTs to attract tourists, constituting motivations for tourism	No. of photographs dominated by NNTs	No. of photographs dominated by NTs	NNT cover in the region	NT cover in the region
	Nature routes from Wikiloc	Geo‐referenced nature routes shared with the public were a mean of assessing society preferences for recreation and tourism	No. of route photographs dominated by NNTs	No. of route photographs dominated by NTs	NNT cover in the region	NT cover in the region
Cultural heritage	Official lists of monumental trees	Monumental trees represent symbols of culture and history, relating to human ‘sense of place’	No. of NNTs in the list	No. of NTs in the list	NNT cover in the region	NT cover in the region
Inspiration	Collective websites of artistic nature photographs	Artistic photographs reflect the choice of inspiring motifs from nature	No. of photographs dominated by NNTs	No. of photographs dominated by NTs	NNT cover in the region	NT cover in the region
Scientific interest	Scopus database of peer‐reviewed scientific literature	The number of scientific publications on NNT or NT species in a country indicates the scientific interest on these species groups in that country	No. of publications on NNTs	No. of publications on NTs	No. of NNT species present in the country	No. of NT species present in the country
Pollen allergenicity	Pollen allergenicity from the database http://allergome.org	The allergenic potential of a tree reduces the benefit of human–nature interactions	No. of NNT species producing allergenic pollen	No. of NT species producing allergenic pollen	No. of NNT species not producing allergenic pollen	No. of NT species not producing allergenic pollen

Ecosystem disservice (EDS)

The contribution of NNTs to the EDS pollen allergenicity was assessed using the same procedure described for CES. In this case, the source of information was the Allergome website (http://www.allergome.org), which compiles worldwide information on allergenicity of plant species. For each of 62 countries/regions we counted the number of NNT and NT species with and without allergenic pollen, using the lists of tree species present in each country. We then calculated the effect size using logOR , in which the observed value was the frequency of NNT species producing allergenic pollen. The expected value was calculated assuming that both NNTs and NTs have the same chances of being allergenic (see Table 2 and Appendix S4B). LogOR values higher or lower than 0, respectively, indicate a higher or lower contribution from NNT species to the EDS than expected by chance.

Predictors of the variation of NNT effects on ecosystem (dis) services

Based on previous knowledge (Brundu & Richardson, 2016; Castro‐Díez et al., 2014a; Kull et al., 2011; Vilà & Pujadas, 2001) we selected a set of nine predictors to explain the variation of NNT effects on ecosystem (dis)services (Table 3). For RES, we selected one biogeographic predictor (biome), two properties of the vegetation structure (native ecosystem type and NNT stand type), one functional property of NNTs (N‐fixing), and the phylogenetic relatedness between species (see Section II.2b). Biome was selected because the literature suggests that the effects of non‐native plants on the nutrient cycle depend on the large‐scale climatic conditions, as captured by biomes (Castro‐Díez et al., 2014a). The type of native ecosystem (e.g. grassland, shrubland, forest, etc., see categories in Table 3) was included to account for the gross functional distance between the NNT and the dominant control vegetation, a key factor to explain the magnitude of the impacts caused by non‐native plants (Castro‐Díez et al., 2014a; Chapin et al., 1996; D'Antonio & Corbin, 2003). The type of NNT stand (whether planted or naturalized) may explain variations in the effects of NNTs on RES because of different functioning between an artificial (planted) stand and a spontaneous self‐maintained system (forest with naturalized NNTs) (Cruz‐Neto et al., 2018; Paz et al., 2015). Finally, the N‐fixing ability of NNTs was selected because of its well‐known effect on soil properties and processes (Castro‐Díez et al., 2014a; Liao et al., 2008; Vilà et al., 2011). For PES, CES, and EDS we also selected biome, plus two indicators of socio‐economic development (gross domestic product and human development index), two demographic predictors (population density and proportion of rural population), and an index of human disturbance (ecological footprint) (Table 3). Socio‐economics and demography are recognized determinants of people's demands for resources and their perception of cultural assets from non‐native plants (Kueffer, 2017; Kull et al., 2011; Vaz et al., 2018; Vilà & Pujadas, 2001). The numeric predictors were not significantly correlated in our data set (Spearman's rank correlation test: r  < 0.045, P > 0.05). Information on the sources used to obtain these predictors and the rationale for their selection is shown in Appendix S7.

Table 3

Predictors used to explain the variation of non‐native tree (NNT) effect size on ecosystem (dis)services across case studies. The last column indicates the category of ecosystem service to which the predictor was applied (RES, regulating; PES, provisioning; CES, cultural ecosystem services; EDS, ecosystem disservice).

Acronym	Description	Predictor categories	Applied to ecosystem service category
(Bio)Geographic context
1. Biome	Biome of the study site or dominating in the country	Tropical forest* Subtropical forest Subtropical desert Mediterranean Temperate forest Temperate grassland/desert Boreal forest	RES, PES, CES, EDS
Stand and species properties
2. Ecosystem	Native ecosystem type	(Semi)desert Grassland shrubland Open forest Forest Urban	RES
3. Stand type	NNTs in forest plantations or naturalized	NNTs in planted stands NNTs in naturalized stands	RES
4. N‐fixation	NNT is N‐fixing or not	NNTs N‐fixing NNTs not N‐fixing	RES
Socio‐economic development
5. GDP	Nominal gross domestic product (US Dollars)	Numeric data	CES, PES, EDS
6. HDI	Human Development Index (ranking values)	Numeric data	CES, PES, EDS
Demography
7. PopDens	Population density (km‐2)	Numeric data	CES, PES, EDS
8. RurPop	Proportion of rural population (%)	Numeric data	CES, PES, EDS
Human disturbance
9. EFP	Ecological footprint (ranking values)	Numeric data	CES, PES, EDS

The term ‘forest’ is used here in a broad sense, including also savannahs and woodlands.

Data analyses

Computation of grand mean effect size

To assess the contribution of NNTs to the provision of each ecosystem service across RES, PES CES and EDS (with the exception of timber provision, see below), all effect sizes (Hedges' d or logOR ) obtained for each ecosystem (dis) service were combined using a random effects meta‐analysis model (REMA) to provide a grand mean effect size [either d + for numerical data (RES) or logOR for frequency data (non‐timber PES, CES and EDS)], where the weight of each case study was the reciprocal of the case study variance. In a random‐effects model, the variance of each study results from the variability within (i.e. sampling error) and among case studies (i.e. the random component). We calculated the latter using the restricted maximum‐likelihood estimation (REML) for numeric data, and the DerSimonian‐Laird (DL) model for frequency data (Borenstein et al., 2009; Viechtbauer, 2010), using the rma() function implemented in the R package metafor (Viechtbauer, 2010). This function also provides the 95% confidence intervals for each grand mean effect size and a two‐tailed parametric test checking whether the effect size differs from zero. Given the non‐normal distribution of the residuals of some models, we additionally assessed the two‐tailed significance of each grand mean effect size through non‐parametric permutation tests under 1000 iterations using the permutest() function from package metafor (Viechtbauer, 2010).

In our data sets some case studies were derived from the same publication (RES), or refer to the same NNT species (RES), or come from the same country (RES, PES and CES), and thus may be more closely related to each other than to other case studies. To explore if non‐independence affected the results, we additionally assessed the grand mean effect size and its significance using multi‐level meta‐analysis models (MLMAs), each including one source of non‐independence (publication, NNT species or country) as random factor (Nakagawa & Santos, 2012). These models allow for different variations of effect sizes between case studies within the same level of the random factor. MLMA models were fitted using the rma.mv() function in metafor. Differences in fit between the REMA and MLMA models were assessed using the Akaike Information Criterion corrected for small sizes (AIC ), so that the more complex models (MLMAs) are considered to be an improvement on the simplest one (REMA) if they result in a reduction of AIC of two points or more (Burnham & Anderson, 2004; Senior et al., 2016).

In the case of timber provision, the effect size of each case study was the difference between the proportions of timber provided by NNTs and by NTs, with no associated variance. Case studies were in most cases the largest administrative regions below the country level (equivalent to European NUTS‐2), but for six countries (Bulgaria, Czech Republic, Chile, Ireland, New Zealand, and Portugal) data were only available at the country level. The grand mean effect size (±95% confidence intervals) was calculated as the weighted median (due to non‐normal distribution) across all case studies. The weighting factor was the proportion of all global timber produced by each country/region (multiplied by 1000 for scaling reasons), so that country/regions with larger annual timber harvest (NNT+NT) contributed more to the grand mean effect size. A two‐tailed Wilcoxon rank test was computed to assess whether the grand mean effect differed from zero, using the wilcox.test() function in R.

Computation of heterogeneity and structured meta‐analyses

For each grand mean effect size calculated with REMA (either d or logOR ), we computed the heterogeneity across effect sizes using several statistics: (i) the Q statistic is the sum of squares of the deviations of each effect size from the grand effect size, weighed by the inverse of the effect sizes' variances. Q was tested against a chi‐squared distribution with n–1 degrees of freedom (n = number of case studies) to assess whether the observed heterogeneity is greater than expected by chance (Borenstein et al., 2009). The caveat of the Q statistic is that its reliability and significance depend on the number of case studies (Borenstein et al., 2009; Nakagawa & Santos, 2012). Thus, we computed two additional statistics; (ii) T which is the estimate of the between‐case study variance; and (iii) I which is the proportion (in %) of the total variation in effect sizes that is due to the between‐study variance (T ) (Borenstein et al., 2009; Nakagawa et al., 2017). The three statistics were calculated using the rma() function of the metafor R package. In the case of the MLMA models, we partitioned I between the random factor level and the case study level following Nakagawa & Santos (2012).

We further assessed whether the variation of effect sizes could be explained by the predictors (fixed factors) shown in Table 3. We performed random‐effects (REMA) structured meta‐analyses using the rma() function, which allows incorporating predictors and returns coefficients and an omnibus test assessing whether the coefficient differs from zero. For continuous predictors, the function also provides the regression slope and its significance. Given the non‐normal distribution of residuals in many cases, we additionally assessed the two‐tailed significance of the predictors over 1000 iterations with the permutest() function (Viechtbauer, 2010).

In the case of RES, we also tested for a phylogenetic signal on the NNT effects using the Phytools R package (Revell, 2012). We first constructed a phylogeny of NNT species, starting from the time‐calibrated molecular phylogeny of Zanne et al. (2014). We selected the taxa of our study (at genus level) using the congeneric.merge() function; then we pruned the phylogeny with the drop.tip() function to obtain a separated NNT phylogeny for each RES with a minimum of 10 species (otherwise the statistical power was considered too low). The phylogenetic signal was assessed using two common comparative metrics (Blomberg's K and Pagel's λ) (Blomberg, Garland & Ives, 2003; Pagel, 1999). Although the rationale behind these metrics is different, both approach 1 when the species phylogenetic signal approximates predictions under a Brownian motion model of evolution, and approach 0 when there is no phylogenetic signal.

To explain the variation of NNT effects on timber provision, separate linear models were conducted with the R function lm() to assess the effects of the biome and socio‐economic predictors in Table 3. In some cases, effect sizes were aggregated at the country level before performing the linear models to match the scale at which predictors were available. Country was tested as an additional predictor, as this data set did not allow us to perform a MLMA with country as a random factor (see Section II.2a). In the analysis of country, the six countries where timber provision was only available at country level were removed from the analyses. To improve homoscedasticity, the effect size (the difference between the proportions of timber provided by NNTs and by NTs) was transformed according to the formula: arcsine (sign(x) × sqrt(abs(x))). To compensate for the different contribution of each case study to global timber production, the weighting factor was included in all linear models with the R function offset().

Publication bias

Meta‐analysis results may be affected by publication bias, i.e. the selective publication of articles finding significant effects over those which find non‐significant effects (Begg, 1994). Publication bias for each RES was investigated by exploring asymmetry in a funnel plot, with effect sizes on the x‐axis and standard error of effect sizes on the y‐axis. In the absence of publication bias, this plot is expected to be a symmetric funnel shape, with a larger dispersion of effect sizes for studies with smaller sample size, i.e. those with large standard errors of effect size (Borenstein et al., 2009). We assessed funnel asymmetry using the random/mixed‐effects version of the Egger's test, which performs a structured meta‐analysis with the standard error as predictor, and returns its slope and significance (Sterne & Egger, 2005). This test was implemented using the regtest() function of the metafor package (Viechtbauer, 2010). A significant relationship implies asymmetry in the funnel plot, which may be an indication of publication bias due to missing values on one side of the funnel. However, there are other reasons for funnel asymmetry besides publication bias, such as heterogeneity (Nakagawa & Santos, 2012; Viechtbauer, 2010). Thus, when the Egger's test on the meta‐analysis without predictors indicated asymmetry, we repeated the test on the meta‐analysis with the predictor which explained more heterogeneity. If this test still reported asymmetry, we assessed the impact of publication bias by removing case studies responsible for funnel asymmetry (Borenstein et al., 2009), and by applying the trim‐and‐fill method (Duval & Tweedie, 2000). This method uses an iterative procedure to remove the most extreme small studies from the asymmetric side of the funnel plot, then adds the original studies back into the analysis, imputes a mirror image for each one, and re‐computes the meta‐analysis. If the new grand mean effect size retains the same sign and significance, then we conclude that publication bias has a trivial or modest impact, but if there is a shift of the sign or significance of the grand mean effect size, then the impact of publication bias may be substantial (Nakagawa & Santos, 2012).

Synergies and trade‐offs between ecosystem (dis)services

Some variables are simultaneously involved in multiple services, which may lead to synergies and trade‐offs among services within and among categories (Cord et al., 2017; de Groot, Wilson & Boumans, 2002). For instance, depth of the litter layer is positively related to soil formation but negatively to fire‐risk prevention (see Table 1), which may lead to trade‐offs between these services. We checked for the existence of positive and negative relationships among the effects of NNTs on different ecosystem (dis)services. Given that the unit of observation differed across types of ecosystem (dis)services (NNT species for RES versus country or administrative region for PES, CES and EDS), associations were evaluated through two separate analyses. For RES, we computed the mean effect size (d ) of each NNT species on each service and tested for significant correlations among the effects of the same NNT species on different RES. For PES, CES and EDS, we calculated the mean effect size (logOR ) for each country/region, and analysed the correlations among the NNT effects found in the same countries/regions on different PES, CES and EDS. We used pairwise Spearman rank correlation tests in R software.

R codes and part of the data used in this study are available in the repository Consorcio Madroño (doi: 10.21950/EGM8SE).

RESULTS

Description of the data set

We analysed a total of 1683 case studies (529 on RES, 154 on PES, 938 on CES, and 62 on the EDS pollen allergenicity). Data covered 44 countries (33 on RES, 22 on PES, 13 on CES and 13 on EDS), all continents except Antarctica, and seven biomes (Fig. 1). Continents with developed countries, such as Europe, were over‐represented with respect to their size, whereas large continents such as Asia and Africa were under‐represented (Fig. S1). Temperate, mediterranean and tropical biomes accounted for most case studies, while other biomes were less (boreal forest, subtropical forest and desert, and temperate grassland/desert) or not (polar) represented (Fig. S2).

Figure 1

Simplified representation of the distribution of case studies. Data were collected to evaluate worldwide effects of non‐native tree species on regulating (RES), provisioning (PES) and cultural (CES) ecosystem services and ecosystem disservices (EDS). For simplicity only RES are represented at the local scale (dots), whereas data for PES, CES and EDS are represented at the country scale (flags). The map shows the biomes considered in this study for illustrative purposes (simplified from the FAO Global Ecological Zones). The term ‘forest’ is used in a broad sense, including also savannahs and woodlands.

The best represented ecosystem services were soil formation (RES), soil fertility (RES), timber provision (PES), recreation and ecotourism (CES), and aesthetics (CES), with more than 115 case studies each. The least represented ecosystem services were pollution regulation (RES), pollination (RES), non‐timber provision (PES) and scientific interest (CES), with less than 20 case studies each (Fig. 2). For plague control and flood regulation (RES), no valid case studies were found.

Figure 2

Effects of non‐native tree (NNT) species on ecosystem services assessed using the random‐effects model (REMA). The mean effect size of NNTs and 95% confidence intervals are depicted across the set of case studies considered for each regulating (A), provisioning (B) and cultural (C) ecosystem services (sample sizes are indicated next to each service). Positive or negative mean effect sizes, respectively, indicate that NNTs (or sites dominated by NNTs) had greater or smaller scores for the service, compared to native tree (NT) species or to control sites dominated by native vegetation. Asterisks to the right of the bars indicate that the mean effect size differs significantly from zero according to a permutation test with 1000 iterations. Values on the right axis indicate the heterogeneity I , which is the proportion (in %) of the total variation in effect sizes that is due to between‐study variance.

Our RES data set covered a total of 125 NNT species. Among them, four species (Ailanthus altissima (Mill.) Swingle, Pinus radiata D. Don, Robinia pseudoacacia L., and Falcataria moluccana (Miq.) Barneby & J.W. Grimes) were the best represented, with 20–22 case studies; four additional species (Acacia saligna (Labill.) Wendl., A. longifolia (Andrews) Willd., Eucalyptus globulus Labill., and Ligustrum lucidum W.T. Aiton) were represented by 10–17 case studies. At the other extreme, 44 NNT species were each represented by a single case study (Table S1). The PES, CES and EDS data sets covered the (dis)services provided by all NNT versus all NT species of a country or region (see Section II.1); thus, the list of NNT species was not specified in this case.

Effects of non‐native tree species on regulating ecosystem services

Our random‐effects meta‐analysis (REMA) of 529 case studies obtained from the scientific literature revealed that NNTs increased climate regulation, soil fertility, soil formation and soil erosion control, but decreased fire‐risk prevention (Fig. 2A). Asymmetry in funnel plots (Fig. S3) and Egger's tests suggest the presence of publication bias for climate regulation, soil fertility and fire‐risk prevention; however, the asymmetry of the latter disappeared when ecosystem type was included as predictor (Table S2). For climate regulation and soil fertility, the trim‐and‐fill procedure did not change the mean effect size, suggesting that the impact of publication bias was trivial (Table S2). In addition, the removal of the six case studies responsible for the funnel asymmetry in the soil fertility data set (see Fig. S3) resulted in a smaller, but still significant, mean effect size (d shifted from 0.63, P = 0.012 to 0.46, P = 0.027). For pollution regulation, pollination and water regulation we found no significant effects of NNTs (Fig. 2A).

Compared with REMA, MLMA including either reference, country or NNT species as a random factor, generally improved the model's explanatory power (AIC reduction ≥2). However, results were largely consistent between REMA and the MLMA. The only exceptions were soil fertility and soil formation, which were not significantly affected by NNTs when NNT species (soil fertility) or country (both RES) were added as random factor in MLMA (Table S3).

The heterogeneity of NNT effects across case studies was high (I  > 80%) for all RES except pollination, Fig. 2A, Table S3). Biome explained the heterogeneity of effects on climate regulation, soil formation and soil erosion control (Table S4), with larger effects in low‐latitude biomes (tropical and subtropical forests) than in middle‐latitude biomes (mediterranean and temperate forest) (Fig. 3A–C). Fire‐risk prevention was more decreased by NNTs in shrublands than in grasslands and forests (Fig. 3D). The increase of soil fertility and formation by NNTs was greater in stands of naturalized NNTs than in NNT plantations (Fig. 3E, F). Finally, higher increases of soil fertility, soil formation, soil erosion control, and water regulation were found when NNTs were N‐fixing (Fig. 3G–J). The remaining effects of predictors on RES were non‐significant (Table S4). The ability of the Fabaceae clade (encompassing N‐fixing species) to contribute most to soil fertility and soil formation was also reflected in the significant values of Pagel's λ, which indicated a phylogenetic signal in these effects. However, for the rest of RES, neither Pagel's λ nor Blomberg's K show a phylogenetic signal attributable to a Brownian model of evolution (Table 4).

Figure 3

Predictors explaining the effects of non‐native tree (NNT) species on regulating ecosystem services (RES) under random‐effects structured meta‐analysis: biome (A–C), native ecosystem type (D), stand type (E, F) and N‐fixation of the NNT (G–J). The figure shows the mean effect size (d ) of NNTs and 95% confidence intervals across the set of case studies considered for each predictor category. Positive or negative mean effect sizes, respectively, indicate that sites dominated by NNTs had a greater or smaller score of the RES than control sites with native vegetation.

Table 4

Results from two common metrics used in comparative analyses (Blomberg's K and Pagel's λ) to test for a significant phylogenetic signal in the effects of non‐native trees (NNTs) on regulating ecosystem services (RES). Each cell contains the value of the metric and its significance (P) according to the expectation of a Brownian model of evolution. N represents the number of NNT species in each RES. Significant results (P < 0.05) are indicated with asterisks.

Metric	Climate regulation	Fire‐risk prevention	Pollution regulation	Soil erosion control	Soil fertility	Soil formation	Water regulation
Blomberg's K (P)	0.025 (0.859)	0.184 (0.289)	0.456 (0.588)	0.056 (0.686)	0.543 (0.091)	0.063 (0.687)	0.087 (0.661)
Pagel's λ (P)	0.000 (1.000)	0.001 (1.000)	0.001 (1.000)	0.117 (0.326)	1.006** (0.001)	0.393* (0.037)	0.001 (1.000)
N	54	35	14	37	56	79	57

Effects of non‐native tree species on provisioning ecosystem services

From our global data set of 144 case studies obtained from worldwide forestry reports, we found that the proportion of timber obtained from NTs was slightly higher than that obtained from NNTs (Fig. 2B). ANOVA tests showed large differences in NNT effects across biomes (F = 6.07, P < 0.001), and across countries (F = 17.91, P < 0.001). Timber production in the subtropical forest biome relied mostly on NNTs while temperate and boreal biomes relied more on NTs (Fig. 4A). Some countries obtained timber mostly from NNTs (Argentina, Chile, Ireland, New Zealand, South Africa, UK), while others relied almost exclusively on NT (Austria, Bulgaria, Canada, Czech Republic, Germany, Japan, Switzerland, USA – with the exception of Hawaii) (Fig. 4B).

Figure 4

Effect size of non‐native trees (NNTs) on timber provision across biomes (A) and across countries (B). For each biome/country the horizontal band represents the median; box limits are defined by the 25th and 75th percentiles; upper whiskers are the smallest of the maximum country/biome value and 75th percentile + 1.5 × box extension; lower whiskers are the largest of the smallest biome/country value and 25th percentile – 1.5 × box extension. Circles indicate extreme values outside the whisker interval. The number of case studies in each biome/country is indicated (biomes/countries with less than three case studies were not included in the statistical analysis).

Among the five uncorrelated socio‐economic and demographic predictors, only the proportion of rural population was negatively related to the effect size (slope = –0.04, P = 0.002), indicating that regions with a higher proportion of rural population rely less on NNTs for timber provision (Table S5).

On the basis of the information collected for 16 countries, we found less uniqueness of non‐timber products in NNTs than NTs, i.e. there are more categories of non‐timber products that are exclusively obtained from NTs than from NNTs (Fig. 2B).

Effects of non‐native tree species on cultural ecosystem services

Our global data set revealed that NNTs have a wide range of effects on the five CES categories considered. According to REMA results, NNTs increased aesthetic values determined from catalogues of plant dealers and inventories of urban parks. For recreation and ecotourism, NNTs were present more often in ecotourism websites, but less in nature routes than NTs than expected by chance (Fig. 2C). NNTs increased cultural heritage (i.e. they were over‐represented in catalogues of monumental trees) but have less scientific interest (i.e. they are the subject of fewer scientific publications, according to the Scopus database) than NTs (Fig. 2C). We found no effects of NNTs on inspiration (assessed from the frequency of occurrence in artistic photographs) (Fig. 2C). Including country as a random factor in MLMA improved all models (i.e. explained a high proportion of the residual between‐case studies variance (I ) and decreased the AIC by at least two units, Table S6), but results were consistent with those of REMA for most CES. The only exception was recreation and ecotourism, where the negative selection of NNTs for nature routes was not significant when country was included as random factor in MLMA (Table S6).

Heterogeneity of effect sizes across case studies was high for most CES (I  > 80%, except for recreation/ecotourism and inspiration, Fig. 2C, Table S6). Biome contributed to explaining the variation of NNT effects on aesthetics and recreation and ecotourism (Table S7). NNTs contributed more to aesthetics in tropical and temperate biomes than in mediterranean and boreal ones, and NNTs were used more in tourism websites in the mediterranean than in other biomes (Table S8). In countries with higher gross domestic products, NNTs contributed more to aesthetics, and were more used in tourism websites, but they were selected less often for nature routes. In countries with higher values of the human development index, NNTs had less effect on aesthetics and recreation and ecotourism services, suggesting larger effects in less‐developed countries. In countries with higher population density or with a higher proportion of rural population, NNTs contributed less to aesthetics but more to recreation and ecotourism. The contribution of NNTs to cultural heritage declined at higher population density. NNTs were less selected as recreation and ecotourism, and inspiration assets in countries with larger ecological footprints (i.e. human disturbance, Table S7).

Effects of non‐native tree species on the ecosystem disservice

Overall, NNTs producing allergenic pollen were not more frequent than expected, either using REMA or MLMA with country as a random factor (logOR  ± SE = 0.093 ± 0.197 and 0.587 ± 0.365, with REMA and MLMA, respectively, P > 0.05, N = 62, Table S6). Nevertheless, NNT effects on pollen allergenicity varied widely with context (I  = 91%). In tropical and temperate biomes, NNTs contributed more to pollen allergenicity than in mediterranean and boreal ones (Table S8). Higher NNT contribution to pollen allergenicity was also associated with countries with higher gross domestic products. Finally, pollen allergenicity was increased less by NNTs in countries with a higher proportion of rural population (Table S7).

Relationships between ecosystem services

The effects of NNTs on several RES were correlated with each other. Most Spearman correlations were positive [i.e. among soil fertility, soil formation and erosion control, as well as climate regulation, and water regulation (R  = 0.41–0.68, P = 0.010–<0.001)]. However, fire‐risk prevention was negatively correlated with other RES (R  = –0.69 and –0.67, P < 0.001 for soil erosion control and water regulation, respectively, Table S9).

The overall effects of NNTs on PES, CES and EDS within countries/regions were somewhat correlated. The most significant was a positive correlation between the two sources of information used for aesthetics, i.e. catalogues of plant dealers and urban parks (R  = 0.89, P > 0.001). Aesthetics was also positively correlated with pollen allergenicity (R  = 0.56–0.65, P = 0.01), suggesting that many ornamental NNTs produce allergenic pollen. Uniqueness of non‐timber products was negatively related to pollen allergenicity (R  = –0.71, P = 0.03). Although marginally significant, countries where the contribution of NNTs to pollen allergenicity was higher tended to show higher scientific interest in NNTs (R  = 0.58, P = 0.06, Table S10).

DISCUSSION

The data set

Our study is the first comprehensive global analysis of the effects of NNTs on ecosystem (dis)services that brings together information from a broad spectrum of scientific subjects, as well as from online sources. Despite the global scope, information was more abundant for developed countries, mediterranean and temperate forest biomes, and for tree species with high societal interest. This geographical and taxonomic bias is well documented in the ecological literature (Hulme et al., 2013; Martin, Blossey & Ellis, 2012; Pyšek et al., 2008; Wilson et al., 2007). Human population and forest cover are unevenly distributed across biomes, and scientific institutions are more abundant in wealthy regions (Wilson et al., 2007). Our RES data set covered one third of the 430 NNT species known to be invasive (Rejmánek & Richardson, 2013), of which just 22 species accounted for half of the case studies. This fact highlights that scientists repeatedly target the few species with known large impacts on social‐ecological systems (Hulme et al., 2013; Pyšek et al., 2008). The coverage of ecosystem services was limited by the need to find measures suitable for pairwise comparisons and by the need to harmonize data across multiple countries. Despite these limitations, this study offers the most complete and up‐to‐date analysis of current knowledge on NNTs worldwide, and allowed us to identify key gaps for future research.

Effects of non‐native trees on ecosystem services

Overall, we found more increases than decreases in ecosystem services attributable to NNTs. This result is consistent with many NNTs having been deliberately introduced to create or enhance particular ecosystem services (Brundu & Richardson, 2016; Evans, 2009; Potgieter et al., 2017). The increase in several RES due to the presence of NNTs may be attributed to their high productivity, a trait often selected for NNT introduction (Richardson, 1998; Woziwoda, Kopec & Witkowski, 2014). Increased productivity of forests may promote climate regulation (via carbon uptake), and soil formation, fertility and erosion control (through higher root growth and/or the supply of more organic matter to the soil) (Evans, 2009; Mori, Lertzman & Gustafsson, 2017). The decrease in both timber and non‐timber PES by NNT was unexpected, given that the purpose of many tree introductions is the supply of particular products (Brundu & Richardson, 2016; Evans, 2009; MacDicken, 2015). In the case of timber, this is because many temperate regions rely almost exclusively on NTs. The lower uniqueness of non‐timber products found for NNTs suggests that many of these species were introduced to increase the quality or quantity of products already supplied by NTs, rather than to produce novel products (Krumm & Vítková, 2016). Regarding CES, the ornamental value of non‐native plants was previously reported to be associated with the human preference for novelty and unusual features, e.g. the colourful or large flowers of Jacaranda and Magnolia, the crown shape and size of Sequoia and Ficus (Kueffer & Kull, 2017; van Wilgen & Richardson, 2012). This may explain the selection of NNTs as aesthetic, tourist and heritage assets (Vaz et al., 2018). However, this does not explain the reduced touristic value of nature routes caused by NNTs. Ornamental NNTs are usually confined to gardens or urban parks. Outside these areas, NNTs can occur because they were planted for non‐aesthetic purposes (e.g. timber production or land reclamation), or as naturalized escapes from plantations. The traits that promote non‐aesthetic uses of NNTs (e.g. fast growth; Richardson, 1998), or those associated with naturalization success (e.g. profuse seed production or resprouting capacity; Castro‐Díez et al., 2011; Richardson & Rejmánek, 2004) can lead to monotonous and homogenized landscapes which are apparently less attractive to users of nature routes. Despite the growing number of studies on non‐native species (Hulme et al., 2013), the scientific interest in NNTs was overwhelmingly lower than for NTs. This may be because ecological interest in non‐native species and their impacts is recent compared to the long history of research on pristine habitats and native species (Hulme et al., 2013; Vaz et al., 2017a).

Given that biodiversity often promotes multiple ecosystem services in forests (Gamfeldt et al., 2013; Mori et al., 2017; Poorter et al., 2015), the increase of ecosystem services by NNTs found here seems at odds with the low biodiversity usually found in NNT‐dominated systems (Gaertner et al., 2011; Pyšek et al., 2012; Vilà et al., 2011). Our results may have been exacerbated by the fact that the baseline for comparison may include degraded or non‐forest ecosystems. Moreover, replacing NTs with NNTs may maximize a particular ecosystem service at the expense of reducing the ecosystem's capacity to provide multiple services (Dickie et al., 2014a; Evans, 2009; van Wilgen & Richardson, 2012), and most studies covered herein focus on single, rather than multiple, ecosystem services. We also note that maximizing some particular ecosystem services may not always be beneficial for society or ecosystem functioning. For example, in the naturally infertile soils of some parts of Hawaii and the South African fynbos, the increase of soil nitrogen driven by the introduction of N‐fixing NNTs is disrupting ecosystem functioning and altering several ecosystem services, such as soil fertility and water supply (Gaertner et al., 2011; Le Maitre et al., 1996; Vitousek & Walker, 1989). Thus, focussing environmental policies on ecosystem services may overlook the intrinsic value of nature and leave biodiversity under‐protected (Dee et al., 2017; McCauley, 2006; Silvertown, 2015).

Variability of non‐native tree effects on ecosystem services worldwide

NNT effects in most ecosystem (dis)services showed medium to high heterogeneity, with I values above 80% in most cases. This high heterogeneity is typical of ecological studies, where effect sizes derive from different species and different contexts (Senior et al., 2016). Nevertheless, we were able to explain part of this heterogeneity. Nitrogen‐fixing NNTs had especially strong effects on soil fertility, in agreement with previous studies (Castro‐Díez et al., 2014a; Liao et al., 2008; Vilà et al., 2011). Compared with non‐N‐fixing NNTs, N‐fixing NNTs may attain a higher production in infertile or degraded soils, explaining their contribution to soil formation (by supplying more organic matter to the soil), and to erosion control and water regulation (by a denser protective cover of the soil and a greater net of roots belowground). Accordingly, the low proportion of N‐fixing NNTs among planted stands (9.1–14.1%), compared with naturalized stands (41.5–51.2%), may explain the larger effect size of NNTs on soil fertility and formation found in naturalized stands. Certain RES were more affected by NNTs in low‐latitude biomes than in temperate biomes, possibly related to the stronger effect of non‐native plants on the nitrogen cycle in benign climates reported previously (Castro‐Díez et al., 2014a). Many NNTs are selected primarily for their high potential productivity (Richardson, 1998; Woziwoda et al., 2014), a trait that underpins many ecosystem services. However, that high potential productivity would be realized more in environments where the favourable period for plant growth is longer (e.g. tropical forests) (Castro‐Díez et al., 2014a).

The degree of phylogenetic relatedness among species showed a small role in explaining variation of NNT effects on RES. This means that closely and distantly related species have similar effects. The only phylogenetic signal was found for the effects of NNTs on soil fertility and soil formation, and is associated with the phylogenetic proximity among N‐fixing NNTs, which had greater effects on these RES.

The great cross‐country and cross‐biome variation observed in the contribution of NNTs to timber provision might be explained by several non‐exclusive arguments. First, regions with high availability or variety of NT species (e.g. USA, Brazil) have less need to introduce NNTs to supply their timber needs. Second, regions with a non‐profitable pool of NT species (e.g. because of prevailing slow growth rates, as occurs in the Mediterranean Basin) may have favoured plantations of profitable NNTs. Third, some former European colonies (e.g. Argentina, Australia, Chile, New Zealand or South Africa) have a stronger tendency to plant NNTs, due to the colonial ethos of ‘national development’, to cultural links with the colonists' home countries or with other colonies, or to the loss of native cultures that were more dependent on native species (Carruthers et al., 2011; Speziale et al., 2012). Thus, a complex interaction of environmental, social and historical factors seems to have shaped species selection for PES (Kueffer, 2017). This also explains the poor role of socio‐economic and demographic predictors when considered alone. Nonetheless, regions with a small proportion of rural population (or more urbanised regions) seem to rely more on NNTs for timber. This suggests that rural societies tend to rely more on native assets, because the populace has a closer connection with the environment or less need of resources (Carruthers et al., 2011; Shackleton et al., 2007; Speziale et al., 2012).

Our results showed that NNTs contributed more to CES in countries with greater nominal gross domestic products. Wealthy regions foster the trade and maintenance of non‐native plants (Gavier‐Pizarro et al., 2010; Humair et al., 2015; Vilà & Pujadas, 2001), and thus their contribution to CES (Vaz et al., 2018). Nevertheless, NNT contribution to aesthetics, and to recreation and ecotourism decreased in regions with higher human development index (i.e. higher life expectancy, education level, and income), a trend previously observed in the Iberian Peninsula (Vaz et al., 2018). This suggests a higher awareness of the risks associated with NNTs, and thus a higher preference for NTs as ornamental and tourism assets in more developed regions (Nuñez & Pauchard, 2010; Vaz et al., 2018). Finally, we found a lower contribution of NNTs to tourism and inspiration services in regions with larger ecological footprint, suggesting that a higher pressure on natural resources makes people prefer NTs as inspirational and tourism assets.

Some predictors that explain the heterogeneity of NNT effects relate to intrinsic properties of individual NNTs (e.g. N‐fixation ability) or to the environments in which they occur (e.g. biome). Other predictors relate more to the social dimension, such as socio‐economy, demography, and cultural background. Whereas ecological conditions are expected to be more or less stable in the long term, the social dimension (e.g. human perceptions, norms and values, social memory, institutions and rules) is more prone to changes in a few generations. This has important consequences for the way in which NNTs affect ecosystem services through time (Kueffer & Kull, 2017; Kull et al., 2018; Shackleton et al., 2016). There are many examples of people embracing introduced species in their practices and traditions in preference to native species used previously (i.e. the ‘shifting baseline syndrome’; Kueffer & Kull, 2017; Nuñez & Simberloff, 2005; Speziale et al., 2012). This poses an additional risk to focusing environmental policies solely on ecosystem services, whose values for people are changeable over short time frames (Silvertown, 2015).

Despite the patterns of variation in NNT effects on ecosystem services revealed here, a large proportion of variation remains unexplained. This suggests that other predictors that were not considered in our analysis are important, such as NNT functional traits (e.g. plant size, seed mass or leaf habit; Castro‐Díez et al., 2014a; Pyšek et al., 2012), or the historical factors accounting for the cultural and trade relations between distant regions through time (Speziale et al., 2012). Other potential predictors may operate at local scales, and thus could not be included in our global‐scale analysis (e.g. silvicultural practices, individual choices, attitudes, and behaviours; Grove et al., 2006; Kull et al., 2018).

Synergies and trade‐offs among the effects of non‐native trees on ecosystem services

The correlations found among most RES can be explained by the fact that some tree traits simultaneously contribute to different RES. In particular, the potential growth rate of trees may underpin many of the identified links across the effects of NNTs on RES: fast‐growing trees can be important carbon sinks and thus contribute to climate regulation, while simultaneously promoting erosion control and soil formation through rapid development of a protective soil cover, enhancing soil organic matter, and contributing to regulation of the water cycle. Yet, trade‐offs across RES may also arise when the same trait contributes positively to some services and negatively to others (Potgieter et al., 2017). Thus, some fast‐growing trees can also increase fire risk by supplying a high quantity of fuel to the system, and because they invest less resources in protection against disturbances (Herms & Mattson, 1992). Regarding CES, we identified a strong correlation between the two sources of aesthetic information (i.e. catalogues of ornamental plant dealers and tree inventories of urban parks) indicating strong consistency between them. By contrast, the lack of correlation between the two indicators of recreation and ecotourism (i.e. official tourism websites and nature routes) indicates that they capture different aspects of ecotourism attraction: tourism websites may tend to highlight the ‘unusual’ or the ‘spectacular’ to attract visitors (e.g. a plantation of sequoias in northern Spain or conifers along the Garden Route in South Africa), whereas users of nature routes appreciate more ‘pristine’ nature dominated by native species (Vaz et al., 2018). We also found a synergy between aesthetics and pollen allergenicity, suggesting that NNTs with aesthetic value may also exhibit traits that promote pollen allergenicity (e.g. wind pollination, which was found to increase with urbanization; Williams, Hahs & Vesk, 2015). This result converges with others showing high allergenicity in non‐native air‐borne pollen from ornamental NNTs (Belmonte & Vilà, 2004; Bosch‐Cano et al., 2011). Unfortunately, differences in scale of study used for RES compared with PES, CES and EDS precluded us from exploring other synergies and trade‐offs across different categories of ecosystem services.

Methodological limitations and future perspectives

Publication bias in the scientific literature may have affected our RES analysis. Funnel plot asymmetry suggests that the positive effects of NNTs on climate regulation and soil fertility may have been inflated by publication bias. However, the two methods applied (trim and fill and removal of case studies from the extreme of the funnel) suggest consistent results. Future updates might minimize publication bias, e.g. by extending the literature search to grey literature (Borenstein et al., 2009).

Another methodological issue in meta‐analysis is non‐independence among case studies (Nakagawa & Santos, 2012; Noble et al., 2017). Although partly removed by aggregating related case studies, non‐independence remained in our data sets due to multiple case studies derived from the same publication, the same NNT species, or the same country. The impact of non‐independence on results may be assessed by comparing the results from the two meta‐analysis models, MLMA and REMA (with and without the source of non‐independence as a random factor). Only three ecosystem services (soil fertility, soil formation and touristic value of nature routes) changed from being affected by NNTs (REMA) to a non‐significant effect (MLMA with country as random factor, and also with NNT species in the case of soil fertility). This indicates that results for these services cannot be extended to any random set of case studies. For instance, the increase in soil fertility by NNTs was due to N‐fixing trees; thus, a data set with fewer N‐fixing NNTs would be likely to return a non‐significant result. Also, the effect of NNTs on the touristic value of nature routes varied across countries (e.g. it was decreased in Australia and Spain, but increased in Italy and South Africa). Thus, the overall decrease may be attributed to the high proportion of case studies from Australia and Spain (Appendix S6). Future studies should account for this bias by representing each country with a number of case studies proportional to its population or size.

In the case of case of CES, EDS and non‐timber PES, we used a novel indicator‐based approach, previously described in Vaz et al. (2018). This approach allows assessing whether NNTs are preferred or rejected for particular services, using the ‘offer’ (i.e. availability) of NNTs in the region as a reference value. This approach has the advantage of integrating multiple sources of information, allowing reproducibility and updates as the sources expand (Vaz et al., 2018). In addition, sources of information from social media (e.g. tourism, image‐sharing, or commercial websites) have a much wider coverage than traditional scientific media (Richards & Friess, 2015; Wood et al., 2013). Yet, the approach has some caveats: lower data quality (e.g. potential mistakes in species names in catalogues); low resolution of the analysis (effect sizes are calculated for a group of NNTs of a country/region, rather than for particular species); sensitivity of the effect size metric (logOR ) to the choice of data types and control data (Vaz et al., 2018); or the inability to compare magnitudes of logOR across ecosystem (dis)services with data sources provided in different units (e.g. species counts, species cover, or number of photographs). In summary, this novel approach is useful to obtain preliminary insights on the directions of effects of NNTs on certain ecosystem (dis)services assessed poorly by scientific media, and to cover spatial and temporal scales not attainable through traditional scientific methodologies. Future studies should validate our results for specific contexts using other methodologies (Hernandez‐Morcillo et al., 2013).

Our understanding of how NNTs, and non‐native species in general, influence several ecosystem services simultaneously is poor (Vilà & Hulme, 2017). The approach used here to identify associations among responses of ecosystem services to NNTs represents a first crude exploration. Unfortunately, differences in scale of study used for RES (species), and for PES, CES and EDS (countries/regions), together with the limited number of species or countries common to many ecosystem services, precluded us from using a multivariate approach (Spake et al., 2017). Such an approach would allow identifying bundles of services according to their response to NNTs (Raudsepp‐Hearne, Peterson & Bennett, 2010), knowing how the impact on such bundles changes in space (Spake et al., 2017) and thus identifying areas most at risk from NNTs.

CONCLUSIONS

Our comprehensive worldwide review revealed more increases than decreases in ecosystem services attributable to NNTs. However, the strong context dependency of the effects, and the limitations of the concept of ‘ecosystem services’, means that these results must be interpreted with caution. The anthropocentric view of nature prominent in recent conservation literature may encourage the maximization of a few ecosystem services in the short term at the expense of long‐term sustainability and multifunctionality.

Part of the variation in NNT effects across case studies is explained by stable factors (e.g. biogeography, traits of NNTs), but some is due to changeable socio‐economic and demographic factors.

Trade‐offs and synergies between ecosystem services emerge because they may be associated with the same plant traits.

This review revealed some key knowledge gaps that need attention, such as the lack of information from large areas of Africa and Asia, and for many NNT species whose contribution to ecosystem services has not yet been documented.

AUTHOR CONTRIBUTIONS

P.C.‐D. and O.G. produced the study idea; P.C.‐D., A.S.V., O.G., J.S.S., M.V.L., A.A. and A.S. designed the research; all authors contributed to data collection; P.C.‐D., A.S.V. and O.G. organized and analysed the data; P.C.‐D. led the writing and A.S.V., O.G., J.S.S., M.V.L., D.M.R. and M.V. contributed to the first draft; all authors contributed to the final version of the manuscript.

Appendix S1. Sources of information used in each country to obtain the list of native and non‐native tree species.

Appendix S2. List of references identified by our search criteria for regulating ecosystem services.

Appendix S3. How we dealt with pseudo‐replicates in regulating ecosystem services.

Appendix S4. Detailed protocols for statistical analysis.

Appendix S5. Description of data sources of provisioning ecosystem services.

Appendix S6. Description of data sources of cultural ecosystem services and the disservice pollen allergenicity.

Appendix S7. Description of data sources of predictors used to explain the variation of the effects of non‐native trees on ecosystem services.

Fig. S1. Number of case studies on provisioning (PES), cultural (CES), regulating services (RES) and ecosystem disservice (EDS) per (sub)continent.

Fig. S2. Number of case studies on provisioning (PES), cultural (CES), regulating services (RES) and ecosystem disservice (EDS) per biome.

Table S1. Non‐native tree species in the regulating ecosystem services (RES) data set.

Fig. S3. Funnel plots for each meta‐analysis on the effect of non‐native trees on regulating ecosystem services (RES).

Table S2. Analysis of publication bias in the meta‐analyses exploring the effects of non‐native tree species on regulating ecosystem services (RES).

Table S3. Comparison of meta‐analysis models constructed to assess the grand mean effect sizes of non‐native tree (NNT) species on regulating ecosystem services (RES).

Table S4. Heterogeneity (Q ) of effect sizes of non‐native tree (NNT) species on regulating ecosystem services (RES) across case studies explained by four qualitative predictors under random‐effects structured meta‐analysis.

Table S5. Results of the linear models relating the effect size of non‐native tree species on timber provision with five socio‐economic and demographic predictors.

Table S6. Comparison of meta‐analysis models constructed to assess the grand mean effect size (and its 95% confidence intervals, CI) of non‐native tree (NNT) species on cultural ecosystem services (CES) and on the disservice (EDS) pollen allergenicity.

Table S7. Structured meta‐analysis assessing the contribution of predictors to explaining the heterogeneity of non‐native tree (NNT) species effect size on cultural ecosystem services and on the disservice pollen allergenicity.

Table S8. Effect size of non‐native tree species on cultural ecosystem services (CES) and on one disservice (pollen allergenicity), separated by biome.

Table S9. Pairwise Spearman correlation coefficients among the effects of non‐native tree species on different regulating ecosystem services.

Table S10. Pairwise Spearman correlation coefficients among the effects of non‐native tree species on provisioning and cultural ecosystem services, and on the disservice pollen allergenicity.

Click here for additional data file.