Functional Divergence Drives Invasibility of Plant Communities at the Edges of a Resource Availability Gradient

: Invasive Alien Species (IAS) are a serious threat to biodiversity, severely affecting natural habitats and species assemblages. However, no consistent empirical evidence emerged on which functional traits or trait combination may foster community invasibility. Novel insights on the functional features promoting community invasibility may arise from the use of mechanistic traits, like those associated with drought resistance, which have been seldom included in trait-based studies. Here, we tested for the functional strategies of native and invasive assemblage (i.e., environmental filtering hypothesis vs. niche divergence), and we assessed how the functional space determined by native species could influence community invasibility at the edges of a resource availability gradient. Our results showed that invasive species pools need to have a certain degree of differentiation in order to persist in highly invaded communities, suggesting that functional niche divergence may foster community invasibility. In addition, resident native communities more susceptible to invasion are those which, on average, have higher resource acquisition capacity, and lower drought resistance coupled with an apparently reduced water-use efficiency. We advocate the use of a mechanistic perspective in future research to comprehensively understand invasion dynamics, providing also new insights on the factors underlying community invasibility in different ecosystems.


Introduction
Biological invasions are one of the most important causes of biodiversity loss [1,2]. Currently, there are no signals suggesting that this process is going to be hampered [3]; since alien species have new opportunities to establish and spread outside their native ranges fostered by both socioeconomic factors [4,5] and climate change [6]. Invasive Alien Species (hereafter IAS) severely affect natural habitats and species assemblages [7][8][9], and have also been claimed as one of the main factors driving ecosystems towards biotic homogenization [10][11][12].
The functional features promoting species invasiveness and community invasibility (sensu [13]) have been frequently analyzed in the last decade [14,15], but no consistent empirical evidence emerged so far on which traits or trait combination may promote habitat or community invasibility [16,17]. Most studies in this field have focused on comparative analysis of functional traits between native and invasive species (e.g., [18]), showing that IAS are characterized by a 'fast-return strategy' [19], that is, higher values of functional traits associated with performance (i.e., photosynthetic rate, growth rate, size) and lower values of cost-associated traits (i.e., specific leaf area, wood density) [18], which convey higher competitive ability leading to invasion success. In addition, IAS also showed smaller trait ranges and larger trait divergences compared to native species [20]. However, several studies invoked that a multiple suite of traits may promote IAS spread in different environments [21,22], suggesting that the invasion process is strongly context-dependent [23].
If invasion success is driven by large-scale environmental effects ('environmental filtering'), it might be hypothesized that most successful invaders share similar functional adaptations as resident native species [24][25][26]. Alternatively, if there is a great functional difference between native and invasive species pools, it may be inferred that IAS occupy a different trait niche than co-occurring species, in order to avoid niche overlap and competition with resident species by using untapped resource pools ('niche divergence'; [27][28][29][30]). Additionally, it is not obvious that a species with greater invasiveness becomes a successful invader, as community invasibility seems to be strongly linked to the abiotic and biotic features of the recipient communities [30,31]. Other key issues to consider when dealing with invasion success are phylogeny, intraspecific trait variation, the spatial scale of investigation, and resource availability [17,32]. The latter plays an important role in the invasion process, since high-resource ecosystems tend to be more invaded thanks to greater availability of light, water, and soil nutrients, among others; even though high levels of invasion occur also in some low-resource ecosystems [33].
Recently, functional diversity intended as the variation of traits among organisms [34] has gained a pivotal role in assessing the impacts of alien species on resident native communities [35]. This facet of biodiversity reflects the functional role of species within a community, and its alteration may threaten ecosystem functioning and stability. Novel insights on the functional features promoting the vulnerability of communities to biological invasion may arise from the analysis of 'mechanistic functional traits' (sensu [36]), like those associated with plant hydraulics and water relations. These traits have been seldom included in trait-based studies, probably because of their difficult and time-consuming measurement, but their effectiveness in explaining broad patterns of evolution and ecology in plants has already been proved [37]. Mechanistic traits, unlike general or soft functional traits (which often represent 'syndromes' that can be driven by different physiological functions, e.g., Specific Leaf Area), are clearly linked to distinct physiological functions and can provide deeper insight on species' physiological features, being more strongly associated with habitat characteristics and plant functions with respect to trait syndromes [37][38][39][40].
To the best of our knowledge, only a few studies aimed at disentangling functional diversity patterns using mechanistic traits (see, for instance, [41] for a subtropical forest; [42] for a coastal ecosystem), but these have never been used to test for functional niche overlap between native and IAS and, in particular, to infer how the functional space determined by native species could influence community invasibility in two sites characterized by a different degree of resource availability. In this study, we aimed at (i) determining the functional strategies of native and invasive assemblage (i.e., environmental filtering hypothesis vs. niche divergence) and (ii) assessing the combination of community-weighted functional traits of native assemblage related with vulnerability to the invasion of natural plant communities.

Study Area and Sampling Design
Two sampling sites were selected at the extremes of a resource availability gradient in northeastern Italy ( Figure S1). The first sampling site (Site 1) is the 'Doberdò lake' (centroid coordinates 45.83111 °N-13.56173 °E; total area ~4 km 2 ), a wetland dominated by species adapted to mesic and flooded environments such as Populus nigra and Salix spp. along with wet meadows dominated by Carex spp. Wetlands are often characterized by higher levels of nutrients especially if surrounded by rural areas due to runoff or enriched groundwater, making them prone to biological invasions ( [43] and references therein). The second sampling site (Site 2) is represented by the coastal site of Marano and Grado Lagoon (centroid coordinates 45.72216 °N-13.24836 °E, total area ~163 km 2 ), characterized by sand dunes and saltmarshes habitats. Coastal environments such as sand dunes are characterized by a steep ecological gradient in a relatively restricted space due to the presence of harsh abiotic conditions such as marine aerosol, sand blast, and lower levels of nutrients which strongly shape local plant communities [44]. In addition, coastal areas are experiencing high levels of invasion across the world [45,46]. Both sites are subjected to anthropogenic disturbances (e.g., tourism, proximity to urban areas), even though greater effects are observable in the coastal site. Due to the presence of species and habitat of high conservation value, both sites belong to a network of European natural protected areas (Natura 2000 Network). Climate in the study area is sub-Mediterranean, with rainy cool winters and long and relatively dry summers. Mean annual temperature averages 15 °C, and mean annual precipitation sums up to about 1000 mm (source https://www.meteo.fvg.it, reference period 1999-2018 accessed on February 19 th , 2019). Due to the different morphology and characteristics of the study sites, we collected plant species abundance estimated with visual cover in each sampling unit using a slightly different approach in the two sites but maintaining the same grain size (i.e., squared units of 16 m 2 ). In Site 1, data were collected by means of a stratified random sampling where the number of sampling units was displaced proportionally to habitat area (n = 22). In Site 2 (n = 131), which is characterized by a complex mosaic of small habitat patches and a steep ecological gradient, we used random belt transects since their effectiveness in sampling coastal systems has already been proved [47]. Transects were displaced with variable length (from 16 to 168 m) according to dune extension and coast morphology (see Tordoni et al. [42] for more details). Vascular plants occurring within each plot were identified following Pignatti [48], and nomenclature was standardized with the Taxonomic Name Resolution Service ( [49]; http://tnrs.iplantcollaborative.org/). Furthermore, plant species were classified according to their invasion status using the updated checklist of the Italian alien vascular flora [50]; all alien species have been included in the analysis irrespective of their invasion status.

Measuring Functional and Mechanistic Traits
Functional traits were measured on at least one individual for each species within each site. Leaves or entire individuals were collected, wrapped in cling film, put in humid sealed plastic bags, and stored in cool bags until processing in the laboratory. Fieldwork was carried out in springsummer 2017/2018. Functional traits measured in this study reflect the trade-offs between carbon investment for leaf construction and potential photosynthetic carbon gain and are typically associated with the 'leaf economics spectrum' (LES, [51]). Mechanistic traits related to the efficiency of water transport within leaves (leaf venation architecture), water-use efficiency (leaf isotopic composition), and drought resistance (leaf water potential at turgor loss point) were also measured. Specifically, the following leaf traits were considered: leaf mass per area (LMA, mg cm −2 ), leaf dry matter content (LDMC, mg g −1 ), minor vein length per unit area (VLAmin, mm mm −2 ), water potential at turgor loss point (Ψtlp, MPa), N and C content (N, %; C, %), C to N ratio (C:N), and C and N stable isotope composition (δ 13 C, ‰; δ 15 N, ‰).
LMA and LDMC are considered as proxies for plant resource use efficiency [52]. Specifically, LMA estimates the leaf-level cost of light interception [53], with lower values of LMA usually associated with higher photosynthetic rates and higher nutrient contents [54]. LMA and LDMC were calculated as follows: Fresh leaves were first rehydrated overnight, and leaf turgid weight was measured with an analytical balance along with their area using the software ImageJ [55]. Leaves were then oven-dried for 48 h at 70 °C before measuring dry weight.
Leaf venation architecture comprehends several structural features influencing plant performance, and higher VLA is generally associated with higher leaf hydraulic conductance and gas exchange rates [56]. We measured the length per unit area of minor veins (VLAmin) as: To measure VLAmin, fresh leaves were cleared in NaOH 1M for 48-72 h at room temperature, carefully replacing solution when it turned from transparent to dark-colored. After initial clearance, small portions of leaves of about 1 cm 2 were cut and bleached in NaClO 5% for 1-2 min. Then, samples were treated in a sequence of ethanol solutions at increasing concentrations (25%, 50%, 75%, 100%) and maintained in an alcoholic solution of toluidine blue (3%) overnight. Finally, samples were processed in a series of ethanol solutions at decreasing concentrations and microscopic slides were prepared. Images of small portions of leaves (~5 mm 2 ) were captured with an optical microscope (4× magnification) equipped with a digital camera (model Syrio-2, Pbinternational, Mumbai, India) and VLAmin was measured using PhenoVein software [57].
Ψtlp is strongly associated with species' drought tolerance, and more negative values are generally associated with a greater capability of maintaining turgor under water shortage [58]. Whole individuals were rehydrated overnight by emerging roots in tap water. Then, one or few leaves from each individual were excised and roughly crumbled before being sealed in cling film and immersed in liquid nitrogen (LN2) for 2 min. Leaves (still sealed in cling film) were carefully ground and stored in sealed plastic bottles at −20 °C until measurements. Samples were thawed at room temperature for 5 min before measuring π0 with a dew point potentiometer (Model WP4, Decagon Devices Inc., Pullman, Washington, USA). Ψtlp was then calculated according to [59].
δ 13 C is a proxy for water-use efficiency (WUE), whereas δ 15 N is related to soil nitrogen availability as well as to plant resource acquisition and use [60,61]. Similarly, elemental composition in terms of N %, C %, and C:N provides information on carbon and nitrogen investment costs which, in turn, are correlated with leaf photosynthetic capacity [62]. Leaves were oven-dried (70 °C for 48 h) and then pulverized in a mortar. Dried and ground samples were analyzed for carbon and nitrogen contents (% dry weight) and carbon and nitrogen stable isotope composition by means of an elemental analyzer/continuous flow isotope ratio mass spectrometry using a CHNOS Elemental Analyzer (vario ISOTOPE cube, Elementar, Hanau, Germany) coupled with an IsoPrime 100 mass spectrometer (Isoprime Ltd., Cheadle, UK). All isotope and elemental composition analyses were performed by the Center for Stable Isotope Biogeochemistry (University of California, Berkeley). Long-term external precision based on reference material 'NIST SMR 1577b' (bovine liver) is 0.10‰ and 0.15‰, respectively, for C and N isotope analyses.
The whole list of functional traits included in this study along with their functional significance is summarized in Table 1.

Functional Characterization of the Plant Communities
Firstly, plant community functional compositions for native and alien assemblages were characterized separately at the plot scale using Community-Level Weighted Means (CWMs, [63]) calculated as follows: where wip is the relative abundance of species i in plot p and ti is the mean trait value of species i, computed using the R package 'FD' [64]. In order to assess which functional composition of native assemblage is more related with community invasibility, this was expressed using the normalized alien cover (Inorm) at plot level calculated as: where Calien is the alien cover in the plot and Ctotal was obtained summing up native and alien species cover. This procedure was necessary to exclude eventual empty spaces and make CWM values independent from Inorm [42]. Functional composition niches (i.e., region of the functional space containing all the trait combinations displayed by a selected sampling unit) were computed using the trait probability density approach (TPD, [65]), which reflects the probability of observing different trait values in a given unit (e.g., species, community, region), thus incorporating the probabilistic nature of functional niches. We computed TPDs using CWMs of native and invasive assemblages in each sampling site through the function TPDs in R package 'TPD' [66]; this approach was necessary due to the contemporary lack of intraspecific trait data for most of the species along with the lack of some functional trait measurements for others. We then computed overlap-based dissimilarities that quantify functional differences between these assemblages. In detail, functional dissimilarity (β0) was expressed as 1-Overlap where: where TPDij is the joint density distributions of native and invasive assemblage, respectively. β0 ranged between 0, when two units share the same functional features (complete functional overlap), and 1 otherwise (no functional overlap). To pinpointing the ecological processes relying on TPDs, we decomposed β0 in two complementary components following Carmona et al. [65]. The first reflects the portion of functional volume which is uniquely occupied by one of the units (PU), whereas the other one represents the amount of relative abundance in functional trait space included, or nested, in the joint volume (PN). Community invasibility was analyzed by means of a Linear Mixed-effect Model (LMM) using the 'nlme' package [67]. Inorm was modeled as a function of native CWM for each trait (fixed effects); the site was specified as the random effect in order to control for its variability and for the different sampling design. Furthermore, to account for the effect of spatial autocorrelation in parameter estimation, a Gaussian correlation matrix was added incorporating spatial structures of the data. Starting from a full model including the CWMs of all measured traits, we carried on a backward stepwise variable selection procedure through AICc minimization criteria in order to find a Minimum Adequate Model (MAM) using package 'MuMIn' [68]. R 2 values were calculated using the 'r2glmm' package [69]. This statistic was based on the standardized generalized variance approach, which is the proportion of generalized variance explained by the fixed predictors. Prior to analysis, quantitative CWMs were standardized to have a mean of 0 and unit variance. All statistical analyses were performed using R 3.6.2 [70].

Results
Overall, 73 species were sampled of which 59 natives and 14 aliens (12 invasive species and 2 naturalized). Three IAS were shared between the two sites namely Ambrosia artemisiifolia, Amorpha fruticosa, and Xanthium orientale (for the full list of species, see Table S1). We observed significant differences in terms of native species richness and community invasibility between the two sites, whereas differences in IAS richness were not significant (Figure 1).

Figure 1.
Boxplots of the variation of native, invasive species richness and community invasibility in the two sampling sites. Wilcoxon rank-sum tests were also reported. IAS = Invasive Alien Species.
TPDs of native and invasive assemblages are reported in Figure 2, where a net distinction of functional niches (β0) can be observed between native and invasive assemblages for most of the CWMs in both sites. However, this difference was not always consistent across functional traits and sampling sites. The TPDs of LMA, LDMC, and Ψtlp showed contrasting patterns in terms of the different portions of the functional space occupied by native and invasive species between sites ( Table 2). Across traits, Site 2 showed greater dissimilarities between IAS and native species spanning from 0.29 (C:N) to 0.80 (C); in contrast, β0 ranged between 0.81 (C) and 1.00 (N) in Site 1.
AICc-weighted standardized regression coefficients revealed that the main determinants of community invasibility were generally consistent with the minimum adequate models, having also a comparable effect size (R 2 ≈ 0.42), and suggesting that the functional features promoting invasibility were consistent in both sites. Mechanistic traits of native species related to WUE and drought resistance (δ 13 C and Ψtlp, respectively), and leaf construction costs (C:N) were among the most important factors in determining community invasibility, even though it is worth noting that LDMC emerged as a significant factor in all the selected models ( Table 3). The MAM minimizing the AICc value (best fit 67.54, R 2 = 0.42) showed that Inorm has a positive relationship with LDMC and Ψtlp, whereas a negative relationship was observed with δ 13 C and C:N, even though these factors accounted for a different proportion of variance explained ( Table 4), suggesting that resident native communities more prone to invasion are those with higher resource acquisition capacity, and lower drought resistance coupled with an apparently reduced WUE.

Discussion
To the best of our knowledge, this is one of the first studies investigating possible relationships between community invasibility and mechanistic traits, with a focus on physiological features related to drought resistance and WUE. Among the hypotheses generally advanced to explain community invasibility, namely environmental filtering vs. niche divergence, we provide evidence supporting the latter, even though there is a variation in this signal across sampling sites and functional traits (CWMs; see Figure 2). The most interesting pattern emerging from our data is that invasive species displayed, on average, higher functional niche dissimilarity in both sampling sites, though Site 2 showed lower β0 (see Site 2 in Figure 2, Table 2), in agreement with the hypothesis that native and invasive species pools need to have a certain degree of differentiation in order to coexist in highly invaded communities [71]. Although both sampling sites displayed similar IAS richness, coastal areas were clearly more invaded, in line with previous studies in north-eastern Italy [45,72] as well as on a global scale [46]. The greater community invasibility in Site 2 may be explained by pulses of resources associated with reduced competition ('fluctuating resource theory', [73]) along with constant propagule pressure from adjacent environments subjected to anthropic disturbances, which might act as invasion epicenters. Indeed, a large number of nutrients were carried from inland through river discharges which is usually offset by the water exchange with the adjacent open sea [74]. Moreover, there might be an effect of intraspecific trait variation which could increase niche breadth and benefit interspecies competition [32,75]. Table 2. Functional dissimilarity (β0) and overlap between native and invasive assemblages in the two sampling sites for each functional trait (CWM) used in this study. PN and PU represent the relative components of the dissimilarity, which is exclusive or joint between the two distributions, respectively. See Table 1 for the full list of abbreviations of the measured traits.

Site 1
Site 2 IAS were often reported to be more performant than native species displaying higher values for traits related to resource acquisition such as photosynthetic capacity and specific leaf area [18,42,76], and this signal was detected for most of the traits measured in this study (i.e., higher peaks of TPDs in invasive assemblage in LMA, LDMC, δ 13 C). Greater phenotypic plasticity, the ability to spread over long distances, and functional features that were not yet present in the native resident community are usually considered as key-factors which may promote a greater community invasibility [20,77], even though there is still an open debate since these factors may vary across time, and invasion stages [78]. According to biotic resistance, native communities should repel less competitive invaders (for instance through reduced survival rates) or invaders sharing similar functional features; however, it has been reported that IAS overcome biotic resistance by competitioninduced trait shifts [79]. In addition, fine-scale experimental studies [79,80] indicated that functionally distinct plant species might have a better chance to succeed in the invasion process being better competitors or avoiding direct competition and co-occurrence with closely related natives, occupying different spatio-temporal niches [17,81]. Perhaps, a great sharing of functional niches between the two subcommunities may hamper community invasibility (e.g., through competition), thus newly introduced species will likely become naturalized rather than invasive [77]. Interestingly, different niche dissimilarity values were observed between sites highlighting once again the strong contextdependence of invasion patterns [22,77] and how these mechanisms are species-specific (few IAS are shared between the two sites, see Table S1).
In stress-prone ecosystems such as coastal dunes and saltmarshes, plant communities tolerate strong limiting factors (e.g., sand burial, salinity, and drought; [44,45]), which probably act as filter selecting only those IAS which are able to withstand them. As a consequence, the greater niche overlap observed in such stressful environments (both in terms of total overlap than of PN component) is not completely surprising, as well as the larger range of the CWMs (see for instance C, δ 13 C, and LMA), probably due to the mosaic of different habitats present in this site in a relatively restricted space, which often host rare and endemic species with peculiar functional attributes [82]. In contrast, Site 1 showed a net separation in the use of the functional space (higher β0) between the two assemblages suggesting that greater niche availability may correspond to different exploitation of the functional space and, in turn, of resources. Generally, it is expected to observe a reduction in trait range with increasing environmental severity [83,84], but this pattern may be inverted if species use contrasting strategies to deal with stress [85,86].

Functional Features of Native Species Pool Facilitating Community Invasibility
In this study, we highlighted the suite of functional traits (CWMs) of native species which were associated with greater levels of invasion. First of all, our results showed that the best performing models (i.e., those with reduced AICc) include mechanistic traits (see Table 3), thus advocating their use in future research. This outcome was somehow expected, since it has been already demonstrated that mechanistic traits are strongly associated with distinct physiological functions and habitat characteristics much more than general 'syndrome-like' functional traits [39,40], even though we found a quite unexpected role of LDMC in our models both in terms of frequency and explanatory power (see partial R 2 in Table 4). Table 3. Multi-model inference from Linear Mixed Models (LMMs) Inorm (response variable) as a function of the Community-Weighted Means (CWMs) of the native assemblage and site as a random factor. Model selection was based on values of the Akaike Information Criterion corrected for small sampling sizes (AICc) across models with all possible factor combinations. Please note that only the models with ΔAICc ≤ 2 were reported along with their R 2 statistics. Secondly, we pointed out how community invasibility depends critically on the functional characteristics of the recipient community [87]. Indeed, we observed that communities more susceptible to invasion are those in which the most abundant native species, on average, invest more resources in leaf construction costs, are less resistant to drought, and have a greater stomatal aperture. According to the LES [51], native species invest more biomass to build 'conservative' leaves which make them more resistant to drought stress and promote longer leaf lifespan and higher survival chances under abiotic and biotic stress [76]; in contrast, IAS are characterized by an 'acquisitive' strategy which minimizes leaf construction and maintenance costs [19,33]. Specifically, native communities displaying a 'conservative strategy' (e.g., higher LDMC and lower Ψtlp) are more invasible than communities on the opposite of the LES continuum, as observed also in Catford et al. [78]. Regarding drought resistance, our model showed that community invasibility is greater when water availability is not a limiting factor, and this can influence especially seedling establishment [88]. According to the 'fast-slow economic spectrum' proposed by [89], species with lower values of Ψtlp also have higher water transport efficiency and consequently higher C flux and higher growth rates; on the other hand, they are more vulnerable to extreme events and resource fluctuations and thus more prone to invasion [73]. Furthermore, higher levels of invasion were detected in plots where the native assemblage showed lower values of δ 13 C, highlighting that IAS may also obtain a competitive advantage due to a higher WUE [18]. Nevertheless, in our study, this relationship is surely partially driven by the presence of a few species displaying C4 metabolism in Site 2, which is characterized by intrinsically higher values of δ 13 C [90].

Model logLik AICc
Overall, the most successful IAS are expected to be pre-adapted to the new environment, especially in low-resource/highly selective ecosystems, as highlighted by the larger trait ranges in Site 2. Resource availability may be one of the main constraints influencing the establishment of new invasive species which is usually able to better exploit an increase in resource availability than native species [91]. Biotic interactions (i.e., competition and facilitation) may become more and more important along the subsequent stages of the invasion continuum [31]. One of the main consequences of the invasion processes could be the disruption and homogenization of the native resident community [92][93][94][95][96], resulting in a loss of species and a reduction in the functional space [42]. Probably, accounting for phylogenetic signal and intraspecific trait variation along with the use of long-term experiments [78] may help to comprehensively understand invasion dynamics, even though not all traits are phylogenetically conserved.
Supplementary Materials: The following are available online at www.mdpi.com/1424-2818/12/4/148/s1, Figure  S1: a) Map of the study area showing the location of the two sampling sites; right inset displays their position with respect to Italian peninsula; b) Details of the two sampling sites (right panel Site 1, left panel Site 2); red dots and lines represent the plots and the transects, respectively. Table S1: List of vascular plant species sampled in the study area. *denote invasive alien species, † denotes alien species which are naturalized. In bold the IAS shared between the two environments.
Author Contributions: E.T. and F.P. collected, and analyzed the data. G.B. and A.N. conceived the idea. E.T. led the writing of the manuscript; all authors have read and agreed to the published version of the manuscript. We also thank two anonymous reviewers for providing insightful suggestions to improve the manuscript. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.