Proportion of non-native plants in urban parks correlates with climate, socioeconomic factors and plant traits

Urban Forestry & Urban Greening 63 (2021) 127215 Available online 12 June 2021 1618-8667/© 2021 Published by Elsevier GmbH. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). Proportion of non-native plants in urban parks correlates with climate, socioeconomic factors and plant traits ´ Alvaro Bay´ on a , Oscar Godoy b , No¨ elie Maurel c , Mark van Kleunen c , d , Montserrat Vil` a a , e , * a Estaci´ on Biol´ ogica de Do˜ nana (EBD-CSIC), Av. Am´ erico Vespucio 26, 41092, Sevilla, Spain b Department of Biology, University of Cadiz, Cadiz, Spain c Ecology, Department of Biology, University of Konstanz, Konstanz, Germany d Zhejiang Provincial Key Laboratory of Plant Evolutionary Ecology and Conservation, Taizhou University, Taizhou, 318000, China e Department of Plant Biology and Ecology, University of Seville, Seville, Spain ARTICLE INFO Handling Editor: W Wendy McWilliam Keywords: Hardiness zones Invasive species Ornamental plants Plant trait profile ABSTRACT Urban parks and gardens provide cultural and aesthetic services critical for human well-being. Yet, they represent one of the main reasons for the intentional introduction of ornamental species, some of which can escape and establish in natural ecosystems. Besides aesthetic reasons, climate and socioeconomic factors can also modulate which species are planted in urban parks. Here, we evaluate the relationship between traits of 486 ornamental woody species from 46 Spanish urban parks and climatic and socioeconomic variables. We specifically assessed how plant traits, climatic, and socioeconomic factors are related to the proportion of non-native species and, among them, to the proportion of established non-native species in Spain. Overall, we found clear associations between species traits and climatic variables. Most notably, parks with warmer winters have more plant species with conspicuous flowers, whereas parks with colder winters and a more continental climate have more species with higher tolerances to cold and shade. Most of the species recorded in our study are non-native (82 %). Higher proportions of non-native species in urban parks were positively associated with towns with large size homes and inhabitants with higher median age but negatively related to parks with species with a higher hardiness-zone range. Moreover, a greater proportion of non-native species that can establish in the natural ecosystems was found in parks with lower continentality conditions. Our results show that Spanish urban parks have an overwhelming proportion of non-native woody species, some of which have the potential to establish, and that the variation in their proportions can be explained by climatic, and socioeconomic factors. 1. Introduction Urban parks and gardens provide cultural and aesthetic services critical for human well-being (Bolund and Hunhammar, 1999; Hulme, 2007; Kendal et al., 2012), yet they are one of the main reasons for the intentional introduction of ornamental non-native plant species worldwide (Mayer et al., 2017; van Kleunen et al., 2018). This source of non-native species introduction is partly explained by the fact that in our continuous search of novelty, we, as humans, are actively looking for plant features that increase the aesthetic value of urban parks while creating a sense of place. In addition, previous work has shown that such search is also modulated by the socioeconomic context of the towns within a country. In general, the presence of non-native plant species is explained by climate, human factors and a combination of both (Pyˇ sek et al., 2010). For instance, towns with high average incomes and high development dedicate greater monetary investment to the introduction of novel species with higher aesthetic value (Vaz et al., 2018) and to the maintenance of green areas, with more diverse species composition (Pyˇ sek et al., 2010). Likewise, densely populated regions are characterized by higher propagule pressure, therefore, more individuals are planted and more invasive plants present (Pino et al., 2005; Pyˇ sek et al., 2010). Nevertheless, the decorative value is only one reason for planting non-native species in urban park. Other important reasons relate to the ability of some non-native species to cope with environmental stress. Commonly, park managers select species that are easy to establish and maintain. This means finding species with particular traits that are preadapted to the particular climatic conditions. However, this * Corresponding author. E-mail addresses: [email protected] (´ A. Bay´ on), [email protected] (M. Vil` a). Contents lists available at ScienceDirect Urban Forestry & Urban Greening journal homepage: www.elsevier.com/locate/ufug https://doi.org/10.1016/j.ufug.2021.127215 Received 30 November 2020; Received in revised form 2 June 2021; Accepted 7 June 2021 Urban Forestry & Urban Greening 63 (2021) 127215 2 management can bring environmental problems as it also facilitates their establishment in natural ecosystems in the introduced range (Gonz´ alez-Moreno et al., 2014; Maurel et al., 2016; Pyˇ sek et al., 2010). In general terms, it has been well-documented that the distribution of woody species along climatic gradients is modulated by the particular organs and whole-plant traits that allow species to tolerate environmental stressors such as frost, shade, and drought (Rueda et al., 2017; Zanne et al., 2014). Therefore, when a non-native species is introduced in a novel area, the interaction of climatic conditions with the species traits determine which species will persist (Haeuser et al., 2018, 2017; Dullinger et al., 2017; van der Veken et al., 2008). This interactive process is not different in the case of urban floras, although gardening practices (e.g. watering, drainage systems for excess water, pruning) can overcome or compensate for some of these climatic limitations. The introduction of non-native woody species in urban parks can pose a conservation problem because some of the species can escape, establish populations (i.e. naturalize) in natural ecosystems and become invasive. Several invasive woody species have the potential to disrupt the functioning of natural ecosystems causing environmental impacts, as well as socio-economic and human health problems (Hulme, 2007; Niinemets and Pe˜ nuelas, 2008; Parker et al., 1999). Although the number of naturalized species has increased exponentially during the last decades (Seebens et al., 2017), we still have a poor understanding of the drivers influencing the introduction of non-native plant species into the urban flora compared to the drivers influencing the latter stages of the invasion process. Along the invasion process from transport to spread and invasion (Blackburn et al., 2011; Colautti and MacIsaac, 2004), it has been well documented that species traits, climatic conditions and socioeconomic factors ultimately influence invasion success. Nevertheless, these factors are also likely to influence the probability that a species is introduced somewhere in the first place, and it is worth noting that the species traits determining invasion success might not necessarily be the same as those determining their deliberate introduction, since the latter depends directly on human perception, preferences and uses. With this knowledge at hand, what remains to be explored is how these three main components (i.e. species traits, climate and socioeconomic context) explain variation in the proportion of non-native ornamental species (non-native species number/total species number) in urban parks. Indeed, the relative importance of each factor is likely to vary across a large territory if there is variation in climatic and socioeconomic variables. For instance, severe climatic conditions might limit the importance of socioeconomic variables in determining the proportion of non-native ornamental species due to direct environmental filtering processes (Kraft et al., 2015). However, with a relaxation of stressful conditions due to mild climatic conditions or indirect effects mediated by human gardening activities, human population size or their income is expected to influence the proportion of non-native ornamental species planted in urban parks. Critically, these extrinsic factors are going to be mediated by plant traits (Vaz et al., 2018). Therefore, knowledge of the trait profiles of ornamental plant communities in concert with climatic and socioeconomic correlates will increase our understanding of the processes driving the introduction and composition of non-native floras in urban parks, and the establishment of these non-native species. In these lines, previous work has shown that socioeconomic factors explain as much of the variation of the distribution of non-native species established in natural ecosystems as climatic factors (Essl et al., 2011; Pyˇ sek et al., 2010). However, we do not know to what extent the proportion of woody species naturalized in natural ecosystems is associated with median species trait values of the parks, and their climatic and socioeconomic variables. In this study, we analyze the species planted in urban parks that vary in size and founding year. We specifically ask (1) how climatic and socioeconomic characteristics of the town in which the parks are located, determine the parks’ plant-trait profiles and (2) how plant origin (native vs. non-native) and the invasion status (established vs. not established in natural ecosystems) of the species planted in parks are associated with plant traits, climatic and socioeconomic variables? To answer these questions, we focused our study on the peninsular territory of Spain (hereafter Spain) for several reasons. First, the ornamental use of nonnative plant species in Spain is one of the main pathways of introduction of potential invasive species (Sanz Elorza et al., 2004). With this detailed information, we can assess the trait profiles of the communities of non-native urban floras that are more likely to escape from parks and become invasive in natural landscapes. Recently, 83 invasive or potentially invasive species were identified in nursery catalogues, available for commerce, including eight regulated invasive species (Bay´ on and Vil` a, 2019). The existence of these invasive and potential invasive species for consumers make it easier to be planted in parks and gardens despite that it is forbidden. Second, Spain has a broad range of climatic conditions (from Mediterranean and semiarid to temperate and continental climates). We expect that inland cold conditions will impose a strong filter to non-native woody species, while drought effects in some coastal areas in southern Spain will be ameliorated by irrigation during summer. Moreover, the drought periods characteristic for most climates in Spain has evolutionary limited the trait profile of native flora to develop species with a summer flowering period and dense foliage (Godoy et al., 2009), and we would expect this to be reflected in non-native species. Consequently, we expect a broader variation of trait profiles of non-native plants in coastal towns. Finally, the towns selected are not equally wealthy across the country. In general, the wealth of a town and its unemployment rates might act as important drivers differentiating the non-native floras of urban parks (Schwartz et al., 2006; Trentanovi et al., 2013; Vaz et al., 2017). 2. Methods 2.1. Urban parks dataset We obtained floristic catalogs of 46 urban parks and gardens distributed in 23 towns across the Spanish peninsular territory (Appendix 1) from the ‘Vivirlosparques’ database (www.vivirlosparques. es/; last accessed 12 Feb 2021), an online tool by the ‘Spanish Association of Public Parks and Gardens’ (AEPJP, 2010). This database provides the most updated and complete information of the species planted within urban parks in Spain as well as additionally important details such as total park area, the founding year and the spatial coordinates for most parks. Missing information of park area was completed by using Google Maps tools and, of the year the park was founded, by city hall websites. From these catalogs, we selected the complete list of 486 woody plant species (i.e. trees, shrubs, arboreal cacti and palm trees). Herbs were excluded because their plantings can vary radically from year to year. Once the species list was obtained, we used The Plant List (2013) to harmonize scientific species names (www.theplantlist.org; last accessed 5 October 2018). Non-specific taxa, such as genera with unspecified species epithets, such as Rosa spp. or hybrids such as Citrus ×aurantifolia, were excluded. In the case of infraspecific taxa (e.g. subspecies, varieties), such as Cedrus atlantica var. glauca, we only kept the binomial species names. Then, each species was classified, according to their origin, as native or non-native in Spain, based on Sanz-Elorza et al. (2004) and Castroviejo (2012). For non-native species, we consulted their establishment status in Spain in (Sanz Elorza et al., 2004), which follows the Richardson et al. (2000) status definitions: established and not established in natural ecosystems. We included as established species, those that can be found growing in natural ecosystems in Spain. These include naturalized species that have established persistent populations, but we here also included non-natives classified as casuals (i.e. species that are frequently found in natural areas but do not attain persisting populations). We included casual non-natives among the established ones because the naturalized vs casual classification is not very precise for introduced species in Spain, and because many species ´ A. Bay´ on et al. Urban Forestry & Urban Greening 63 (2021) 127215 3 considered casual are known to cause environmental impacts in natural areas (Andreu et al., 2009). 2.2. Socioeconomic and climatic dataset Climatic variables collected from WorldClim 2 (Fick and Hijmans, 2017) were calculated for each park with QGIS software (QGIS Development Team, 2009). Raster maps with 1 km resolution were obtained for all the variables available in the source: 11 temperature variables and 8 precipitation variables, as well as altitude (Table 1a). Finally, we obtained information of numerous socioeconomic variables for the total area of the 23 towns selected. Specifically, we obtained all the variables available at the Spanish Statistics Office (INE, 2017), which provide information for mean values between 2010 and 2016 (Table 1b) of 11 variables associated with population demography, 4 with urbanization, and 6 with employment and economy. We decided to obtain all the information available because we do not have ‘a priori’ expectation of which variables are relevant for our study. Such high dimensionality was later reduced when performing statistical analyses (see below). Due to the limited resolution of the database of socioeconomic Table 1 Climatic and socioeconomic variables used to describe 46 urban parks in 23 towns across Spain. A) Climatic variables Type of variable Variable Code Kind of variable Unit First PCA selection Temperature Annual mean temperature mean_t quantitative Celsius degree (◦C) * Mean diurnal range diurn_rang quantitative Celsius Isothermality isotherm quantitative % Temperature seasonality t_season quantitative % Maximum temperature of the warmest month max_t_wm quantitative Celsius Minimum temperature of the coldest month min_t_cm quantitative Celsius * Annual range of temperature t_ann_rg quantitative Celsius * Mean temperature of the wettest quarter mt_wett_Q quantitative Celsius Mean temperature of the driest quarter mt_drie_Q quantitative Celsius Mean temperature of the warmest quarter mt_warm_Q quantitative Celsius Mean temperature of the coldest quarter mt_cold_q quantitative Celsius Precipitation Annual precipitation ann_precip quantitative mm Precipitation of the wettest month prec_wettM quantitative mm Precipitation of the driest month prec_drieM quantitative mm Precipitation seasonality prec_seaso quantitative % Precipitation of the wettest quarter prec_wettQ quantitative mm Precipitation of the driest quarter prec_drieQ quantitative mm * Precipitation of the warmest quarter prec_warmQ quantitative mm Precipitation of the coldest quarter prec_coldQ quantitative mm Altitude alts quantitative m a.s.l. B) Socioeconomic variables Type of variable Variable Code Kind of variable Unit First PCA selection Basic demography Number of residents n_resid_10000 quantitative x 10 000 Proportion of residents ≤14 years old n_0.14y quantitative % Proportion of residents 15−64 years old n_15.64y quantitative % Proportion of residents ≥65 years old n_.64y quantitative % Median age median_age quantitative years * Proportion of national residents nationals_percent quantitative % Proportion of foreign residents foreign_percent quantitative % Proportion of residents native from the town native_percent quantitative % * Proportion of residents born outside the town born_outs_percent quantitative % Natality rate natality_per_thousand quantitative ‰ Mortality rate mortality_per_thousand quantitative ‰ Urbanistics Number of homes n_homes_10000 quantitative x 10 000 Number of conventional homes (Catastral) conv_homes_10000 quantitative x 10 000 Mean household size (number of inhabitants per home) home_size quantitative * Proportion of single-person homes single_pers_home_percent quantitative % Employment and economics Unemployment rate unemploy_percent quantitative % Working population in active age (20−64 years old) X20.64y_percent_ocupated_in_active quantitative % * Active population activity_percent quantitative % Industry workers industry_workers_percent quantitative % Service workers service_workers_percent quantitative % Average money incoming per home avg_inc_euro_10000 quantitative x 10 000 € * Area of the town area_ct quantitative km 2 C) Park-dependent variables Type of variable Variable Code Kind of variable Unit Park-dependent Area of the park area_m_sq quantitative m 2 Founding year of the park year quantitative Table 1. A first compilation of climatic variables (A) was performed through a systematic search from WorldClim 2 (Fick and Hijmans, 2017). The first compilation of socioeconomic variables (B) was performed through a systematic search from the Spanish Statistics Office (INE, 2017). Park-dependent variables (C) were mainly obtained from AEPJP (2010), and secondary, from Google Map tools (area) and city hall websites (year). Finally, we made a selection using a Pearson correlation analysis (Appendix 2, Fig. A2.1 and A2.2) and PCA (Fig. 1a and b) in which we applied a broken-stick criterion (Appendix 2, Fig. A2.4 and A2.5), as described in the Methods. * Variables marked with an asterisk were finally selected by the PCA and broken-stick method. ´ A. Bay´ on et al. Urban Forestry & Urban Greening 63 (2021) 127215 4 variables, possible variations below the town level, such as districts or neighborhoods, were not possible to address. Therefore, the socioeconomic data of parks within the same city were the same. We obtained data from a wide range of town sizes —from very densely populated (e.g. Barcelona) to very small towns (e.g. Ciudad Real)— and climatic conditions —from a wet temperate climate (e.g. Bilbao) to a hot and dry Mediterranean climate (Lorca)—. The average number of parks within the same city is 2, and only 2 towns have more than 3 parks (Appendix 1). 2.3. Plant traits dataset We used several sources of information to obtain organ-level and whole-plant traits at the species level. The majority of plant traits were extracted from the TRY Database (Kattge et al., 2011). We performed a second search in several internet databases (Table 2). For each species, we obtained 21 plant traits (17 continuous and 4 binary) related to whole-plant characteristics, leaves, reproductive organs (flowers and dispersal units), plant flowering and fruiting phenology, and tolerances to several environmental stressors such as frost (minimum absolute temperature tolerated and number of frost days per year), shade, drought, waterlogging and fire, as well as the minimum and maximum hardiness values and its range (Table 2). Hardiness zones are cartographical standardizations of the average annual minimum winter temperatures divided into 10-degree F zones (S´ anchez de Lorenzo-C´ aceres, 2004; USDA, 2012). Thus, the hardiness values of a given species is defined as the ability to survive and grow in a specific range of hardiness zones. We consider the lower hardiness zone value that the species can tolerate as “minimum hardiness”, the highest value as “maximum hardiness”, and the different between both values, the “range of hardiness”. Spain includes hardiness zones 7a (-17.8 ◦C to -15 ◦C) to 11b (7.2 ◦C–10 ◦C; Fick and Hijmans, 2017; S´ anchez de Lorenzo-C´ aceres, 2004). This comprehensive search did unfortunately not render trait information for all species. Therefore, we excluded for further analyses five traits with data missing for more than 50 % of the species (Table 2). For the remaining traits, missing data was imputed. Previous work has shown that imputing missing values is better than simply removing missing-data variables or missing-data species, as removing data could bias the results (Penone et al., 2014). We applied nonparametric missing value imputation using random forests implemented in the R package ‘missForest’ (Stekhoven and Buhlmann, 2012), which uses a random forest trained on the observed values to predict the missing values. This imputation method can handle multivariate data consisting of continuous and categorical variables simultaneously, as was the case in our study. Overall, after seven iterations, the out-of-bag normalized root mean squared error was low (1.1608; PFC 0.0295), indicating a good imputation performance (Stekhoven and Buhlmann, 2012). 2.4. Statistical analyses 2.4.1. Overview Due to the high number of variables obtained, and some of them showing high collinearity, we decided to reduce the number of variables using statistical tools with the final objective to avoid overfitting the statistical models (King and Jackson, 1999; Zhou et al., 2010). We first built a database of the plant-trait profile of the parks, calculating the median value of each quantitative trait and the mean value of each binary trait (frequency) for each park. We considered each park as a replicate. The reduction of variables was done on each of the three main groups of variables selected (i.e. climatic variables, socioeconomic variables, and park trait profile) following three complementary steps. First, we performed a correlation matrix to identify the most correlated variables (Appendix 2). If two variables were closely correlated (r-Pearson >0.70), we removed one of them. Second, we performed a principal component analysis (PCA), in which we used a broken stick criterion (Simon et al., 2011) on the first and second dimensions (Appendix 2) to identify the main variables explaining the obtained PCA. Third, we performed Lasso regression models (Friedman et al., 2010; Simon et al., 2011; Tibshirani, 1996) to exclude non-informative explanatory variables (Hesamian and Akbari, 2019). With the completion of these three steps, we obtained the final set of variables, which were used to best predict variation in the proportion of plant origin (native vs. non-native) and in the invasion status (established vs. not established in natural ecosystems). This statistical process is fully explained in Appendix 3. According to this procedure, the first and second axes of the PCA on climatic variables explained 39 % and 35 % of the variation, respectively (Fig. 1a). We selected five climatic variables (annual mean temperature, precipitation in driest quarter, minimum temperature in coldest month, and annual range of temperatures; Table 1a). For the PCA on socioeconomic variables, the first and second principal axes explained 37 % and 25 % of the variation, respectively (Fig. 1b). We selected five socioeconomic variables (mean household size, median age of inhabitants, percent of native population, percent of working population between 20 and 64 years old in active, and average incoming money per home; Table 1b). Finally, the first and second principal axes of the PCA involving park trait profiles explained 28 % and 19 % of the variation, respectively (Fig. 1c). We selected five plant traits (minimum, maximum and range of hardiness, tolerance to shade, and having conspicuous flowers; Table 2). With the final set of variables selected, we then performed generalized linear models (GLMs), using the ‘glm’ function from the R package ‘stats’ (R Core Team, 2019) with a linear combination of previously selected variables as predictors in each case, the species richness of the parks as prior weights, and area and founding year of the park as covariates because bigger and older parks have had more space and time to accumulate a wider diversity of species. 2.4.2. Relationships between plant-trait profile and climatic and socioeconomic characteristics To answer the first question, namely whether climatic and socioeconomic variables explain variation in the functional plant-trait profile across parks, we considered in the Lasso regressions all the climatic and socioeconomic variables preselected as independent variables, while we treated each preselected trait (i.e. minimum, maximum and range of hardiness, tolerance to shade and having conspicuous flowers) as dependent variable in each model (see ‘Overview’). After the reduction of variables, we performed five GLMs (Poisson distribution with the quasipoisson setting to account for overdispersion), one for each trait selected of the park profile. Predictors were a linear combination of climatic and socioeconomic factors previously selected in the Lasso regression (i.e. minimum temperature in the coldest month, annual temperature range, and household size; Table A4.1 from Appendix 4), and area and founding year of the park were included as covariates. 2.4.3. Influence of plant-trait profile, climatic and socioeconomic characteristics on the proportion of non-native species To answer the second question, namely, whether species origin: native vs. non-native; and establishment status: established vs. not established; can be explained by climatic and socioeconomic variables, and the parks’ plant-trait profile, we first obtained the number of native and non-native species within each park, and for the non-native ones, the number of species that have established in natural ecosystems (being invasive, naturalized or casual, according to Sanz Elorza (S´ anchez de Lorenzo-C´ aceres, 2004; USDA, 2012)(2004) and other works (Bay´ on and Vil` a, 2019; Gass´ o et al., 2010). We used a paired Student’s t-test (‘t. test’ function of the ‘stats’ basic R package; R Core Team, 2019) to test whether the numbers of non-native and native species were similar. Similarly, we tested whether the numbers of not established and established non-native species were similar. ´ A. Bay´ on et al. Urban Forestry & Urban Greening 63 (2021) 127215 5 Table 2 Traits used to describe woody plant species in 46 urban parks in 23 towns across Spain. (B¨ artels and Schmidt (2014); Bot´ anica (2014); Brines et al. (1997); Centraal Bureau voor de Statistiek (2003); CSIC (2014); Czerepanov (1995); Gallagher and Leishman (2012); Hintze et al. (2013); Moles et al. (2004); Morales (2002); Plants For A Future (2010); S´ anchez de Lorenzo-C´ aceres (1999); Universidad de M´ alaga (2011); USDA (2008); Wheatley (2009) and Wyse-Jackson (2006)). Type of trait Trait Code Kind of variable Unit % real data % imputated data First PCA selection Data sources beyond TRY Whole plant Max height height_m quantitative m 100 % 0 % 1, 2, 3, 4, 5, 6, 7, 9, 10, 11, 13, 14, 16 Max lifespan — quantitative years 44 % Var.Ex. 1, 2, 3, 4, 5, 6, 7, 11 Leaves Area leaf_area_sq_mm quantitative sq. mm 54 % 46 % 1, 2, 3, 4, 5, 6, 7, 10, 17 Lamina length leaf_length_cm quantitative cm 85 % 15 % 1, 2, 3, 4, 5, 6, 7, 10, 17 Lamina width leaf_width_cm quantitative cm 63 % 37 % 1, 2, 3, 4, 5, 6, 7, 10, 17 Leaf type leaf_type_0to1 binary 0 =broadleaved 97 % 3 % 1, 2, 3, 4, 5, 6, 7, 10, 18 1 =conifer Reproductive organs Floral unit size — quantitative mm 9 % Var.Ex. 1, 2, 3, 4, 5, 6, 7, 10 Conspicuous flower visual_flower_0to1 binary 0 =No 93 % 7 % * 1, 2, 3, 4, 5, 6, 7, 10 1 =Yes Pollination syndrome pollin_synd_0to1 binary 0 = Anemophilous 93 % 7 % 1, 2, 3, 4, 5, 6, 7, 10, 12 1 = Entomophilous Gymnosperm flowers gymnosp_flower_0to1 binary 0 =No 93 % 7 % 1, 2, 3, 4, 5, 6, 7, 12 1 =Yes Dispersal unit mass — quantitative g 15 % Var.Ex. 1, 2, 3, 4, 5, 6, 7, 10, 12, 17 Plant Phenology Reproductive phenology timing phenol_timing quantitative ordinal month 85 % 15 % 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 14, 16 Leaf phenology leaf_phen_0to1 binary 0 =deciduous 97 % 3 % 1, 2, 3, 4, 5, 6, 7, 8, 9, 11, 15, 18 1 =evergreen Tolerances To frost temperature tol_frost_temp quantitative Celsius 82 % 18 % 1, 2, 3, 4, 5, 6, 7, 9 To frost days — quantitative num. days 32 % Var.Ex. 1, 2, 3, 4, 5, 6, 7 To shade tol_shade_0to5 quantitative 1−5 92 % 8 % * 1, 2, 3, 4, 5, 6, 7, 9, 11 To drought tol_drought_0to5 quantitative 1−5 89 % 11 % 1, 2, 3, 4, 5, 6, 7, 9, 11 To waterlogging tol_waterlogging_0to5 quantitative 1−5 76 % 24 % 1, 2, 3, 4, 5, 6, 7, 9 To fire — quantitative 1−5 29 % Var.Ex. 1, 2, 3, 4, 5, 6, 7, 11 Hardiness min. hardiness_min quantitative 1−12 67 % 33 % * 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 Hardiness max. hardiness_max quantitative 1−12 67 % 33 % * 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 Hardiness range hardiness_range quantitative 1−12 Calculated directly: difference between max. and min. values * Calculated Data sources Num. Source reference 1 Brines R, Tejuelo I, Beltr´ an P and Balaguer ´ A. (1997) Guía Verde. https://www.guiaverde.com/. Accessed 1 Sep 2017 2 CSIC Arbolapp. (2014) In: Arbolapp. http://www.arbolapp.es/. Accessed 1 Sep 2017 3 Morales J. (2002) Infojardín. http://www.infojardin.com/. Accessed 1 Sep 2017 4 S´ anchez de Lorenzo-C´ aceres JM. (1999) ´ Arboles Ornamentales. http://arbolesornamentales.es/. Accessed 1 Sep 2017 5 Universidad de M´ alaga. (2011) Jardín Bot´ anico de la Universidad de M´ alaga. http://www.jardinbotanico.uma.es/jardinbotanico/index.php. Accessed 1 Sep 2017 6 Wyse-Jackson P. (2006) Missouri Botanical Garden. https://www.missouribotanicalgarden.org/. Accessed 1 Sep 2017 7 Bot´ anica Y Jardines. (2014) In: Bot´ anica Y Jardines. http://www.botanicayjardines.com/. Accessed 1 Sep 2017 8 B¨ artels A. and Schmidt P. A. (2014) Enzyklop¨ adie der Gartengeh¨ olze, 2nd edn. Ulmer E. Verlag 9 Plants For A Future (PFAF). (2010) In: Plants For A Future. https://pfaf.org/user/Default.aspx. Accessed 1 Jun 2017 10 Weathley, R. (2009) B and T World Seeds. http://b-and-t-world-seeds.com/. Accessed 1 Jun 2017 11 USDA. (2008) USDA Plants Database. https://plants.usda.gov/java/. Accessed 1 Jun 2017 12 Hintze C, Heydel F, Hoppe C, et al (2013) D3: The Dispersal and Diaspore Database – Baseline data and statistics on seed dispersal. Perspectives in Plant Ecology, Evolution and Systematics 15:180–192. 13 Centraal Bureau voor de Statistiek BioBase (2003) In: Centraal Bureau voor de Statistiek. https://www.cbs.nl/nl-nl/onze-diensten/methoden/classificaties/overig/bio base-2003/biobase-2003. Accessed 1 Jun 2017 14 Czerepanov SK (1995)Vascular plants of Russia and adjacent states (the former USSR). Cambridge University Press, New York 15 Zanne AE, Tank DC, Cornwell WK, et al (2004) Three keys to the radiation of angiosperms into freezing environments. Nature 506:89–92. 16 Moles AT, Falster DS, Leishman MR, Westoby M (2012) Small-seeded species produce more seeds per square metre of canopy per year, but not per individual per lifetime. Journal of Ecology 92:384–396. 17 Gallagher RV, Leishman MR (2012) A global analysis of trait variation and evolution in climbing plants. Journal of Biogeography 39:1757–1771. 18 Panchen ZA, Primack RB, Nordt B, et al (Panchen et al., 2014) Leaf out times of temperate woody plants are related to phylogeny, deciduousness, growth habit and wood anatomy. New Phytol 203:1208–1219. Table 2. A first compilation of plant traits was performed through a systematic search in the TRY Database (Kattge et al., 2011). A second search was performed to complete the database —see Data Sources. Then, those variables for which we had more than 50 % of the data were selected, and, on them, we made a nonparametric missing value imputation by random forest (Stekhoven and Buhlmann, 2012). ‘Var.Ex.’ indicates variables excluded from imputation because not enough data. Finally, we made a selection using a Pearson correlation analysis (Appendix 2, Fig. A2.3) and PCA (Fig. 1c) in which we applied a broken-stick criterion (Appendix 2, Fig. A2.6) as described in the Methods. * Variables marked with an asterisk were finally selected by the PCA and broken-stick method. ´ A. Bay´ on et al. Urban Forestry & Urban Greening 63 (2021) 127215 6 Fig. 1. PCAs for relations of variables among parks. Principal component analyses (PCA) for A) climatic variables for each park, B) socioeconomic variables of the localities for each park and C) median values of the plant traits in each park. See Appendix 1for the list of localities and codes of parks, and Tables 1 and 2 for the description of traits and variables. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article). ´ A. Bay´ on et al. Urban Forestry & Urban Greening 63 (2021) 127215 7 Then, for each park we calculated the following proportions of species: •Proportion of non-native species, pAi=Ai Ni+Ai •Proportion of established non-native species, pEi=Ei NEi+Ei Where Ni is the number of native species in park i; Ai is the number of non-native species in park i; NEi is the number of non-native species that have not established in natural ecosystems for park i; and Ei is the number of non-native species that have established in natural ecosystems for park i. We tested whether these proportions varied significantly across parks using a Pearson’s Chi-squared test (‘chisq.test’ function of the R package ‘stats’; R Core Team, 2019). After the reduction of the number of variables by using correlations, PCA and broken stick criteria, we considered for the lasso regressions all the remaining variables and traits preselected as independent variables, while we treated the proportion of non-native species and the proportion of established species as dependent variables in each model (see ‘Overview’). Then we performed two GLMs (binomial distribution), one for each proportion. Predictors were a linear combination of traits, climatic and socioeconomic factors previously selected in the Lasso regression (i.e. annual temperature range, household size, median age of town inhabitants, percent of working population in active age (20−64 years old with job), maximum hardiness and range of hardiness; Table A4.2 from Appendix 4), and area and founding year of the park were included as covariates. 2.5. Data resources The data underpinning the analyses, reported in this paper, are deposited in the Zenodo repository at https://doi.org/10.5281/zenodo.4095422 (Bay´ on et al., 2020) 3. Results Overall, the database has 486 woody species. On average, there were 53 ±5.14 (mean ±se) species per park and the maximum number was 144 species. Many species appear only in one or a few parks. Each species is in 5 ±0.26 (mean ±SE) parks. More precisely, 180 species (37.04 %) appear in only one park each, 68 species (13.99 %) appear in two parks, and 45 species (9.26 %) in three parks, whereas only 11 species (2.25 %) occur in half of the parks or more (Appendix 1). The list of most common species and its status were: Cupressus sempervirens (32 parks; established), Celtis australis (26 parks; native), Robinia pseudoacacia (26 parks; established), Styphnolobium japonicum (25 parks, nonnative), Cedrus atlantica (24 parks; established), Cedrus deodara (24 parks; established), Magnolia grandiflora (24 parks; non-native), Cercis siliquastrum (23 parks; established), Laurus nobilis (23 parks; native), Phoenix canariensis (23 parks; non-native) and Populus alba (23 parks; native). According to the origin, we found 86 native species (17.70 %) and 400 non-native species (82.30 %), 130 (26.75 %) of which have established in natural ecosystems (Table 3). 3.1. Park characterization The correlogram of climatic variables (Fig. A2.1 from Appendix 2) shows that there is a strong correlation between different types of climatic variables. Values of precipitation, except seasonality, showed very high correlations between each other (minimum Pearson’s r =0.49 between precipitation of the driest quarter and the coldest quarter). The first axis of the climatic PCA accounted for 38.5 % of the overall climatic variation across parks, and it was mainly explained by the opposite trend between precipitation and mean temperature in the driest quarter. The second PCA axis accounted for 34.9 % of the observed climatic variation across parks, and it was explained by the same behavior of the minimum temperature in the coldest month and the mean in the wettest quarter (Fig. 1a). Among socioeconomic variables, the highest correlations (Fig. A2.2 from Appendix 2) was observed between the number of residents and number of homes (Pearson’s r =1). Average household sizes in the towns of the parks were directly related to the proportion of kids and adult inhabitants (Pearson’s r =0.62 for kids under 14 years old, and r = 0.60 for adults between 15 and 64 years old), and inversely to the proportion of people over 64 years old (r =-0.70). This proportion of people over 64 is in addition directly related with the median age of the population (Pearson’s r =0.97). Moreover, we found that for the socioeconomic PCA, household size and median age were both distributed along the first axis (36.7 %) but in opposite directions. Median age of the population followed a very similar pattern to the percentage of workers in the service sector within the first socioeconomic PCA axis. Finally, the percentage of employed people, as well as the average income per home, were both distributed along the second PCA axis (25 %) and in the same direction (Fig. 1b). Regarding the parks’ plant-trait profiles (Table 2), some variables were highly correlated (Appendix 2, fig. A2.3). For example, entomophilous pollination was correlated with several other traits such as the presence of conspicuous flowers (Pearson’s r =1), maximum height (Pearson’s r =-0.82) and the tolerated minimum temperature (Pearson’s r =0.52). As expected, there was also a high correlation between cold tolerance and the minimum value of hardiness (Pearson’s r =0.79). Other correlations worth mentioning were those found between the type of leaf (conifer needle – broad leaf) and whether or not species are gymnosperm (Pearson’s r =0.98). Leaf size variables were all closely correlated with each other (Pearson’s r =0.76 for width vs area; 0.76 for area vs length; 0.73 for width vs length). With all this variation, the first axis of the plant trait PCA accounted for 28 % of the overall trait variation and it was mainly explained by an opposite trend between minimum hardiness, hardiness range and the presence of conspicuous flowers. In addition, maximum hardiness and leaf width was also distributed along the first axis of variation but overall had a less statistical weight. Shade and drought tolerance were partly related to the second axis of the PCA in opposite directions (Fig. 1c). 3.2. Relationships between plant-trait profiles and climatic and socioeconomic characteristics After reducing model dimensionality, two climatic variables were retained to explain some variation in the trait profile across parks. These variables were the minimum temperature in the coldest month, and the range of temperatures. The remaining predictors, these are, household size, founding year of the park and park area, did not explain any of the traits selected as dependent variables (Table 4). Specifically, minimum and maximum plant hardiness as well as the presence of conspicuous flowers were positively related to minimum temperature in the coldest month. On the other hand, hardiness range and shade tolerance were negatively related to minimum temperature in the coldest month (Fig. 2a). Similarly, we observed a significant positive relationship between the annual range of temperatures and the minimum and maximum hardiness (i.e. higher maximum hardiness values means plants are less tolerant to cold), while the relationship between the annual range of temperatures and the tolerance to shade was negative. However, there was no significant relationship between the annual Table 3 Number of woody ornamental plant species by status and growth form in 46 urban parks in 23 towns of Spain. Trees Shrubs Palms Arboreal Cacti Total Native 59 26 1 0 86 Non-native not established 154 85 30 1 270 Non-native established 95 30 4 1 130 Total 308 141 35 2 486 ´ A. Bay´ on et al. Urban Forestry & Urban Greening 63 (2021) 127215 8 range of temperatures and the hardiness range, nor with the presence of visual and conspicuous flowers (Fig. 2b). Influence of plant-trait profiles, climatic and socioeconomic characteristics on the proportion of non-native species Overall, there is a greater presence of non-native (39.2 ±4.15; mean ±se) than of native species (14.2 ±1.27; paired t-test; t =7.475; p-value <0.01; Fig. 3a) across parks. For the case of non-native species, the number of species that are not known to have established in natural ecosystems (21.15 ±2.06) is significantly larger than the number of established species (18.02 ±2.27; paired t-test; t =2.4393; p-value = 0.0187; Fig. 3b). The proportion of non-native species and the proportion of established species are both significantly heterogeneous across parks (Chi-squared =132.98, and 148.17, respectively; p <0.01 in both cases). After reducing model dimensionality, only one climatic variable (annual range of temperature), two socioeconomic variables (mean of household size, and median age of inhabitants) and one plant trait of the park profile (range of hardiness) explained the proportion of non-native species. The remaining predictors (percentage of active population, maximum hardiness, founding year of the park, and park area) did not explain any of the proportions selected as dependent variables. The GLM analysis indicates (Table 5) a significant negative relationship between annual range of temperature and the proportion of established species but not with the proportion of non-native species (Fig. 4a). Regarding socioeconomic variables, mean of household size and median age both significantly explained the proportion of nonnative species but not the proportion of established species (Fig. 4b, c). Finally, variation in hardiness range of the park’s profile negatively explained the proportion of non-native species, but it did not explain the proportion of established species (Fig. 4d). 4. Discussion We have a poor understanding of the combination of plant traits, climate variables and socioeconomic factors guiding the planting on non-native woody plant species in urban parks. This knowledge is essential to identify how different drivers influence the success of nonnative species from introduction to invasion. Here, we present insights from a large representation of urban parks across Spain showing that they harbor far more non-native than native plant species. Interestingly, this high proportion of non-native woody species (82.3 %) is not associated with a particular set of taxa. Rather, Spanish parks present a high species turnover, meaning that the vast majority of woody plant species appear in one or very few parks, while very few species are present in the majority of parks. The average percentage of non-native woody species in Spanish urban parks is much higher than in Central European urban floras (Pyˇ sek, 1998), and many other regions (Dangulla et al., 2019; de Freitas et al., 2019; Jha et al., 2019; Templeton et al., 2019), where the average representation of non-native species is between 40 % and 50 %. In line with this finding, the percentage of non-native species available in Spanish nurseries is 76 % approximately (Bay´ on and Vil` a, 2019). The small difference between the percentage in parks and nurseries is probably not significant, and can be due to the fact that the study of Spanish nurseries includes all types of plant growth forms, while this study only focus on woody species. Moreover, many plants in parks are no longer sold in Spanish nurseries (Bay´ on and Vil` a, 2019) such as Buddleja madagascariensis, Casuarina cunninghamiana, Elaeagnus pungens, Fraxinus pennsylvanica, Parkinsonia aculeata, Pittosporum undulatum, Ricinus communis and Schinus molle. We have found that the plant trait composition of urban parks can be explained by a combination of socioeconomic and climatic variables. We observed variation in the proportion of species with conspicuous flowers (which is a direct proxy of the pollination syndrome) across urban parks, and this variation was explained by minimum temperature in the coldest month. Species with entomophilous pollination tend to be poorly tolerant to frost, probably due to the dependence of their pollinating insects on milder temperatures. Meanwhile, wind-pollinated plants have a wider range of tolerance to frost. For instance, oak species are known to be frost tolerant, whereas other species like palms do not show the same adaptations. The latter are predominantly found in towns with oceanic climate, where frost is a rare phenomenon. The pollination Table 4 General Linearized Models on plant traits. A) Minimum hardiness Estimate Std. Error t value Pr(>| t|) (Intercept) 1.5667 0.3541 4.4244 0.0001 *** Minimum temperature of the coldest month 0.0456 0.0071 6.4694 < 0.0001 *** Annual range of temperature 0.0191 0.0048 3.9906 0.0003 *** Mean household size −0.0271 0.0448 −0.6033 0.5497 Founding year of the park −0.0002 0.0002 −1.1550 0.2549 . Area of the park −1.9019E07 1.2108E07 −1.5708 0.1241 B) Maximum hardiness Estimate Std. Error t value Pr(>| t|) (Intercept) 2.1459 0.1665 12.8890 < 0.0001 *** Minimum temperature of the coldest month 0.0180 0.0033 5.4110 < 0.0001 *** Annual range of temperature 0.0082 0.0023 3.6049 < 0.0001 *** Mean household size −0.0216 0.0209 −1.0313 0.3086 Founding year of the park −0.0001 0.0001 −0.9207 0.3628 Area of the park −1.8047E08 5.3817E08 −0.3353 0.7391 C) Range of hardiness Estimate Std. Error t value Pr(>| t|) (Intercept) 1.2575 0.4182 3.0068 0.0045 ** Minimum temperature of the coldest month −0.0288 0.0084 −3.4244 0.0014 ** Annual range of temperature −0.0080 0.0057 −1.3926 0.1714 Mean household size −0.0304 0.0513 −0.5929 0.5566 Founding year of the park 0.0001 0.0002 0.6964 0.4902 Area of the park 6.0555E08 1.2795E07 0.4733 0.6386 D) Tolerance to shade Estimate Std. Error t value Pr(>| t|) (Intercept) 2.0327 0.5283 3.8473 0.0004 *** Minimum temperature of the coldest month −0.0376 0.0107 −3.5268 0.0011 ** Annual range of temperature −0.0320 0.0072 −4.4644 0.0001 *** Mean household size 0.0829 0.0687 1.2063 0.2348 Founding year of the park −0.0003 0.0002 −1.3355 0.1893 Area of the park −1.2227E07 1.6867E07 −0.7249 0.4727 E) Presence of conspicuous flowers Estimate Std. Error t value Pr(>| t|) (Intercept) −3.3205 1.1389 −2.9156 0.0058 ** Minimum temperature of the coldest month 0.0815 0.0220 3.7109 0.0006 *** Annual range of temperature 0.0273 0.0146 1.8666 0.0693 . Mean household size 0.1330 0.1388 0.9586 0.3435 Founding year of the park 0.0008 0.0005 1.5061 0.1399 Area of the park −3.0225E07 3.8541E07 −0.7842 0.4375 Generalized linear models of the five plant traits selected as dependent variables, versus climatic and socioeconomic variables preselected by lasso regression (Table A4.1 from Appendix 4) as independent variables. ** 0.001 <p ≤0.01. *** p ≤0.001, otherwise non significant. Graphical representation in Fig. 2. ´ A. Bay´ on et al. Urban Forestry & Urban Greening 63 (2021) 127215 9 syndrome can change in urban floras depending on their native and nonnative origin. For instance, Aronson et al. (2007) documented that conspicuous flowers, that tends to be beautiful, differ in their pollination requirements; however, we did not find this pattern in our study. Differences between results might be related to the fact that the study of Aronson et al. (2007) was conducted in sites with similar environmental conditions, while our study includes towns with a wide range of climates from continental to Mediterranean, and from inland to oceanic climatic conditions. Besides cold acting as a filter for reproduction syndromes, we found that the tolerance to shade increased with higher minimum temperature in the coldest month and with wider thermal range. In regions with an oceanic climate, such as the northern coast, winters are milder and cloudiness is more frequent and lasts longer, especially in the summer months, than in the rest of Spain. However, in regions where the climate is more continental, as in the plateaus, winters are colder and summers are not only hotter but also very sunny with high evapotranspiration rates (AEMET, 2019). These differences are likely forcing plant species selected to be planted in oceanic climates to be more tolerant to insolation, reducing leaf area to avoid water loss, while plants with larger leaves tend to be distributed in areas with higher water availability and lower irradiance (Holmgren et al., 2012; Markesteijn and Poorter, 2009). Finally, the minimum and maximum hardiness values of the plants were directly related to both minimum temperature in the coldest month and higher annual ranges of temperature. This is indeed an obvious relationship, as hardiness is a proxy of the tolerance to low temperatures (S´ anchez de Lorenzo-C´ aceres, 2004; USDA, 2012), and it represents a primary factor in selecting which species are planted. Specifically, those species with high tolerance to cold, are usually not planted in warm places even if they could survive there. This explanation is supported by the fact that the range of hardiness decreases with the minimum temperature in the coldest month. These results suggest there is a trade-off between species that tolerate shade and those that have conspicuous flowers, as well as with those that have higher hardiness values. We expected that differences across parks between the proportion of Fig. 2. Relationship of five plant traits to (A) the minimum temperature in the coldest month and (B) the annual range of temperature. Generalized linear models (GLMs) of plant traits versus minimum temperature in coldest month and annual range of temperature. The continuous trait variables (i. e. minimum hardiness, maximum hardiness, range of hardiness [scale 1–12], and shade tolerance [scale 1–5]) is shown on the main y axis (left side of the graph), while the scale of the unique binary trait variable (i.e. presence of conspicuous flower), converted to proportion of species having per park, is shown on the secondary y axis (right side). * 0.01 <p ≤0.05, ** 0.001 <p ≤0.01, *** p ≤0.001, otherwise nonsignificant. P-values in Table 4. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article). ´ A. Bay´ on et al.