Exploring the Functional Strategies Adopted by Coastal Plants Along an Ecological Gradient Using Morpho-functional Traits

Coastal dunes are characterised by strong interactions between biotic and abiotic factors along a short gradient from the shoreline to the inland region. We carried out an ecological analysis of the vegetation in a protected area of the Italian coast to evaluate the relationships among species abundance, the occurrence of morphoanatomical traits related to leaves, stems, and roots, and soil variables. Three transects were established perpendicular to the shoreline, with 27 plots distributed in the frontal dunes, backdunes, and temporarily wet dune slacks. An analysis based on community-weighted mean values showed that the pioneer communities of the frontal dunes were dominated by ruderals that are well adapted to the harsh ecological conditions of these environments, showing succulent leaves, high limb thickness values, and low values for leaf dry matter content (LDMC). The backdune vegetation was a mosaic of annual herbaceous and perennial shrub communities showing both ruderal and stress-tolerant strategies (clonality, sclerified leaves, high LDMC values, root phenolics) consistent with less extreme ecological conditions. The dune slack areas were dominated by plants showing adaptations to both arid and flooded environments, such as C4 photosynthesis, amphistomatic leaves, and abundant aerenchyma in the roots. The invasive status, C4 photosynthesis, leaf trichomes, and aerenchyma in the roots were significantly correlated with soil humidity, organic matter content, and pH. These results demonstrate the usefulness of anatomical traits (including root system traits) in understanding the functional strategies adopted by plants. Invasive species tended to occupy plots with high levels of soil moisture, suggesting an avoidance strategy for the harsh environmental conditions of coastal sand dunes. Finally, we suggest including information regarding root systems into coastal monitoring programs because they are directly linked to soil parameters useful in coastal dune management and protection.


Introduction
Coastal sand dunes are unique environments characterised by the close interaction of abiotic and biotic factors that may change dramatically across a few metres, driving a characteristic zonation of plant communities along a gradient from the shoreline to the inland regions (Doing 1985;Sýkora et al. 2004;Isermann 2005;Lane et al. 2008). Salinity, water stress, substrate instability, sand burial, wind abrasion, high temperature, and nutrient limitation are widely recognised as important abiotic factors that impact plant community composition and dynamics of coastal dune ecosystems (Rozema et al. 1985;Hesp 1991;Maun 2009). Mediterranean coastal dunes are characterised by higher water temperature and salinity and more limited tides, waves, and meteorological phenomena compared to those of coastal dunes in environments that are impacted by oceanic storms and hurricanes (Weyl 1970;King 1975). The intrinsic dynamism of coastal dune ecosystems makes them particularly interesting for studies of interactions between abiotic processes and plant functional traits. Analysis of plant functional traits is useful for inferring assembly Communicated by Dennis F. Whigham processes in plant communities, particularly along natural stress gradients, which have been extensively studied in coastal sand dunes (Feagin and Wu 2007;Gallego-Fernandez and Martinez 2011;Chelli et al. 2019). Studies focusing on plant functional traits and European coastal dunes have focused on three aspects. First, functional traits have been examined at the species level to explore plant adaptations to environmental factors (García-Mora et al. 1999;Acosta et al. 2006;Ciccarelli et al. 2009Ciccarelli et al. , 2010Spanò et al. 2013;Retuerto 2013, 2014). Second, others have studied functional traits at the community level by analysing the relationship between functional diversity and the phylogenetic structure of plant communities (Ricotta et al. 2012;Carboni et al. 2013;Brunbjerg et al. 2014). Third, a few studies have analysed coastal dune vegetation with Grime's (2001) CSR (competition, stress-tolerance, and ruderality) scheme (Macedo et al. 2010;Brunbjerg et al. 2015;Ciccarelli 2015). Most of the research conducted on plant functional traits in coastal dune communities has been restricted to the general functional attributes of whole plants (such as life form and growth form, plant height, or clonality) and leaf traits or regenerative attributes (such as dispersal syndrome or seed mass). To the best of our knowledge, there have been no extensive functional analyses of plants in coastal dunes focused on morphological and anatomical traits of stems or roots because analysis of belowground plant structures is timeconsuming and requires specific plant anatomy skills. Because so little is known about the dynamics of belowground plant structures in coastal dunes, we decided to simultaneously explore the relationships among species abundance, trait occurrence, and soil variables. We focused our research on an Italian protected area where coastal dunes are well preserved, and dune zonation is still visible. In this study, through a systematic sampling design based on belt transects from the shoreline to inland, we sampled coastal vegetation and measured both the functional traits of plants and the pedological parameters of soils. We selected numerous morphological traits to explore the general functions of plants, such as growth, space occupancy, competitive ability, photosynthetic rate, ruderality defined as resistance to disturbances (e.g. soil instability, sand burial, wind, and herbivory), and stress tolerance (especially to aridity and high solar radiation, which are typical constraints in coastal ecosystems). Moreover, we also included several anatomical attributes of leaves, stems, and roots, which were useful for characterising coastal plant communities along an environmental gradient in the Atlantic lowland forest (Bona et al. 2020). Determining how functional traits respond to environmental conditions to mediate species distribution and patterns of community diversity is a central question in ecology. Several studies have demonstrated systematic variation in community-weighted mean (CWM) trait values along abiotic gradients (see Muscarella and Uriarte 2015 and references therein). In the present study, we expected to find ruderal species dominating the early successional communities of the shoreline-inland gradient with lower values of specific leaf area (SLA) and leaf dry matter content (LDMC), succulent and hairy leaves, and a welldeveloped root system, and stress-tolerant species in the latesuccessional communities of the backdunes with higher values of plant height, SLA, and LDMC, as well as sclerified leaves (Ciccarelli 2015;Bona et al. 2020). Regarding dune slacks, which are located farthest from the shoreline and are temporarily flooded during autumn and winter, we expected to find plants with adaptations to both flooded and arid environments, such as root aerenchyma and sclerified leaves, respectively. Our aims were (1) to explore how CWM trait values vary along the gradient from the shoreline to inland, confirming our hypotheses; (2) to test whether there are any relationships among species abundance, trait occurrence, and soil variables, and explore the usefulness of using anatomical traits as indicators of plant strategies; and (3) provide a basis for appropriate management strategies for coastal dune ecosystems using analysis of the relationship between plant functional traits and environmental factors.

Study Area
This research was carried out in a protected area called Migliarino -San Rossore -Massaciuccoli Regional Park, which is in the northern part of Tuscany, Italy (Fig. 1). The climate is Mediterranean subhumid (C2), with a mean annual temperature of approximately 15°C and mean rainfall of 800 mm (Rapetti and Vittorini 2012). The study area is part of the Natura 2000 network and includes the site of community importance called the Dune litoranee di Torre del Lago (IT5170001). The plant communities follow a typical shoreline-inland zonation related to an ecological gradient (Ciccarelli 2014(Ciccarelli , 2015Acosta and Ercole 2015). Most of these plant communities correspond to habitat types included in the 92/43/EEC Directive (Habitat Directive), which specifies a list of habitats and species of conservation interest in Europe (European Commission 1992). Frontal dunes, which are located close to the shoreline, are colonised by perennial herbaceous plants that are resistant to sand burial and wind. Behind the frontal dunes, there are interdunal spaces dominated by a mosaic of perennial shrub communities and annual grasslands, mixed with a discontinuous maquis of Juniperus sp. pl. (called stable dunes). Dune slacks are in the farthest part of the shoreline and are flooded, especially during autumn and winter, when precipitation is higher. They are characterised by perennial communities that have different vegetation types depending on the soil moisture, oxygen levels, pH, nutrient content, and annual differences in precipitation (Table S1).

Sampling Design
Three transects, approximately 100 m apart and perpendicular to the shoreline, were established in the study area ( Fig.  1). They ranged from 240 to 270 m in length, according to the morphology of the dune systems. The main advantage of using transects is to capture all the variability among plant communities and soil variables along an environmental gradient; however, a stratified random approach is required to avoid leaving unsampled habitats (Elzinga et al. 2017). Nine 2 m × 2 m plots were positioned along each transect. To quantify plant community composition, three strata (called habitat groups) in each transect, along the shoreline-inland gradient, were identified as follows: (i) "frontal dunes", located between the shoreline and the foredune crest; (ii) "backdunes", located behind the foredune crest in the interdunal grasslands and the stable dunes; and (iii) "slack areas", located in the temporarily humid dune slacks (Table S1). Four plots were placed in both the frontal dunes and backdunes, while one plot was placed in the slack areas depending on the dimension of the area occupied by each habitat group, for a total of 27 plots.

Data Collection
The percentage cover (range 0-100%) of each plant species was visually estimated in each plot as the vertical projection of the aerial part of the plants onto a horizontal surface. Plant nomenclature followed the Portal to the Flora of Italy (http:// dryades.units.it/floritaly/index.php). Soil samples were collected from each plot, and three random subsamples of the first 15 cm of soil were mixed. Each sample consisted of approximately 400 g of soil, which was subjected to the following chemical analyses: quantification of soil moisture content, pH, conductivity, total nitrogen, phosphorous content, and organic matter. All soil analyses followed the protocol described by the SISS (1985). Vegetation samples along with the soil and functional trait data were collected from each of the 27 plots. Fieldwork was conducted between May and June 2019. The dune slacks were dry at the time of sampling.
All species occurring within the plots were considered for the functional trait analysis. Twenty-eight morphological and anatomical attributes of the leaves, stems, and roots were measured (Table 1). Life form and growth form were deduced from the  Barrier to apoplastic inflow of ions (Enstone et al. 2003) literature (Pignatti 2017), but we checked the field. Information about native versus invasive status was retrieved from Galasso et al. (2018). The photosynthetic pathway was derived from online databases such as "eHALOPH" (Santos et al. 2015) and "TRY" (Kattge et al. 2020, and else unpublished data from Günther A. and Wright I. via TRY) and was compared with the anatomical analysis. The clonality and plant height were measured in the field. Other traits were obtained in the laboratory. For species that were present in more than one habitat group, intraspecific trait variability was considered by collecting samples of the given species in each habitat group. The thickness of the epidermis, external periclinal wall of the epidermis, and cuticle was measured only on the adaxial surface of the leaf. If the same species was present in several habitat groups, we collected plant material from that species in each habitat group and considered them different samples. For each species and habitat group, samples were collected from 10 different individuals. Whenever possible, rhizomatous and stoloniferous species were collected far from each other (in the vicinity of the plot) to avoid sampling from the same genet. Adult leaves that were fully expanded and showed no evidence of herbivory between 3°and 6°n odes from the apex were collected. At least 15 leaves per species/habitat group were collected; at least 10 were used for morphological analyses, and three to five were used for anatomical analyses (Pérez-Harguindeguy et al. 2013). Stem samples were collected near the base of the herbaceous species, and branches of approximately 1 cm in diameter were collected from woody plants. Samples of whole roots were collected to obtain the basal region (near the stem) and fine lateral roots (woody plants were sampled from young individuals).
The collection and storage of samples and their processing for the measurement and analysis of the functional traits followed Pérez-Harguindeguy et al. (2013), Ciccarelli et al. (2016), and Klimešová et al. (2019). Plant materials intended for use in anatomical analyses were fixed and stored in 70% ethanol. For the morphological analyses, the materials were moistened, placed in sealed zip-lock bags, and transported in a thermal box to the laboratory, where they were rehydrated in the dark at 4°C overnight (Ciccarelli 2015).
Anatomical sections were made by hand with a disposable razor and subjected to staining or reactions with safranin and astra blue (0.01% safranin and 0.1% astra blue in distilled water; Bukatsch 1972 modified by Kraus and Arduin 1997) to show the cell wall, Sudan III (0.5% in ethanol 80%; Sass 1951) for lipids, Lugol solution (1.5% potassium iodide, 0.3% iodine in distilled water; Johansen 1940) for starch, and ferric chloride (0.3% sodium carbonate, 10% ferric chloride hexahydrate in distilled water; Johansen 1940) for phenolic compounds. The leaf blade was fixed semi-permanently with glycerinated gelatine (7% fish gelatine and 50% glycerine in distilled water). The sections were transversal, in the middle third of the limb (between the median rib and the edge) and at the stem base or branches, the root base, and the absorption zone (when available). Measurements were made under a microscope with an ocular fitted with a calibrated scale using a Zeiss 5 + 100/100 micrometric blade. To measure the cuticle, layers with a positive reaction to Sudan III, including a cutinised wall, were considered.

Data Analysis
We calculated the CWM, following Garnier et al. (2004Garnier et al. ( , 2007. Discrete traits were decomposed into different categories to establish the percentage contribution of each category to the community. Community-weighted mean values were compared using the nonparametric Kruskal-Wallis test with Bonferroni correction for multiple comparisons to verify whether there were significant differences between habitat groups. To explore patterns of plant communities along the seainland gradient, non-metric multidimensional scaling (NMDS) was performed. This technique represents samples in a low-dimensional space by optimising the correspondence between the original dissimilarities and the distances in the ordination (Økland 1996). The species abundance matrix (27 plots per 40 species) was subjected to a square-root transformation to balance the highest and lowest coverage values. A similarity matrix was constructed using the Bray-Curtis similarity index. The Pearson correlation coefficient was also calculated to determine which species were more correlated with the NMDS axes. The same similarity matrix was used to perform a multivariate permutational analysis of variance (PERMANOVA; Anderson 2001) to test the differences in floristic composition of the plant communities at the habitat group and transect levels. PERMANOVA is a routine test for the simultaneous response of one or more variables to one or more factors in an analysis of variance experimental design based on any resemblance measure using permutational methods (Anderson 2001). The variables included in the PERMANOVA design as a source of variation were as follows: (1) habitat group (Ha, fixed with three levels: frontal dunes, backdunes, and slack areas); and (2) transect (Tr, random with three levels: T1, T2, and T3). The significance of each factor in plant species composition was tested using a null model based on 999 permutations. PERMANOVA was also applied to the CWM values.
The contribution of each plant species to the characterisation of plant communities was quantified using similarity percentage analysis (SIMPER; Clarke 1993). SIMPER analyses, based on the Bray-Curtis similarities (Bray and Curtis 1957), were used to identify the species that contributed most to the average between-group dissimilarity for pairs of plots grouped according to the plant community types. All analyses were performed using PRIMER v.7 and PERMANOVA+ (PRIMER-E Ltd, Plymouth, UK; Clarke and Gorley 2015).
We simultaneously analysed the relationships among species abundance, trait occurrence, and environmental variables using RLQ analysis (Dolédec et al. 1996) and fourth-corner analysis (Legendre et al. 1997). First, we checked for multicollinearity among the species traits using the variation inflation factor (VIF) test. Chamaephyte and therophyte life forms, growth forms, native status, clonality, plant height, C3 photosynthetic pathway, sclerenchyma, epistomatic leaves, stomatal protection, stem and root starch, and thickness of the endoderm were eliminated because they exceeded the chosen 10 thresholds, which indicates high multicollinearity between the independent variable and the others, so we excluded them from the analyses. The discrete traits were decomposed into different categories, and finally, we analysed 23 functional traits. The VIF test was also used to check for multicollinearity among environmental variables. In this case, only total nitrogen was >10; therefore, it was excluded from the analyses. Second, three matrices were assembled: species traits (Q (qxm) , 40 species by 23 functional traits), environmental variables (R (nxp) , 27 plots by five soil variables), and species abundance (L (nxq) , 27 plots by abundance of 40 species). RLQ is an extension of co-inertia analysis that searches simultaneously for linear combinations of the variables in Q and linear combinations of the variables in R, maximising the covariance and weighting per L matrix (Dolédec et al. 1996). The preliminary step of RLQ analysis is to analyse each matrix separately. Correspondence analysis (dudi.coa function in R) was applied to the species abundance matrix L. The species trait matrix Q contained both quantitative and categorical variables. In this case, we used the Hill and Smith analysis (dudi.hillsmith function in R), which can perform a multivariate analysis on mixed quantitative variables and factors. For the environmental data, all the variables were quantitative; thus, we applied principal component analysis (dudi.pca function in R). Fourth-corner analysis was used to test the associations between individual traits and environmental variables in a D (mxp) matrix. We computed a permutation test using 999 re-samples to evaluate the significant relationships in matrix D (Dray et al. 2014).
The VIF test, as well as the RLQ and fourth-corner analyses, was performed using R 3.6.1 statistical software (R Development Core Team 2019) with the vegan and ade-4 packages (dudi.coa, dudi.hillsmith, dudi.pca, rlq, and fourthcorner functions).

Plant Communities Along the Shoreline-Inland Gradient
A total of 40 vascular plant taxa distributed in 38 genera and 20 families were recorded in the sampled plots (Table S2). The most represented families were Asteraceae (07) and Poaceae (07), followed by Apiaceae and Brassicaceae (03).
The most common species was Medicago littoralis Rohde ex Loisel., which was found in 63% of the plots, followed by is an alien plant that occurs with self-maintaining populations without direct human intervention, while invasive is an alien plant that occurs with self-maintaining populations without direct human intervention and produces fertile offspring at considerable distances from the parent individuals, thus being able to spread over a large area. Neophytes are alien plants that were introduced in Italy after 1492.
The NMDS bidimensional plot showed that the three habitat groups corresponding to the main plant community types along the shoreline-inland gradient segregated with a stress index of 0.09 (Fig. 2). The first group was composed of the "frontal dune" plots, corresponding to the herbaceous communities in the upper beach and the frontal dunes dominated by Thinopyrum junceum, which was the most important taxon contributing to the total similarity among the plots ( Table 2). The second group, the "backdune" plots, was made up of a mosaic of herbaceous and shrub communities characterised by the dominance of Lomelosia rutifolia and Festuca fasciculata, which contributed almost equally to the total similarity among the plots ( Table 2). The third group, the "slack areas" plots, was characterised by herbaceous communities of the temporarily wet dune slacks dominated by Scirpoides holoschoenus (L.) Soják ( Table 2). The PERMANOVA results showed that the habitat group factor explained a significant (P<0.003) proportion of the variation in species composition at the plot scale (Table 3).
Regarding the CWM of the functional traits (Table 4), the percentage of the contribution of hemicryptophytes was higher in the backdunes than that in the other habitat groups, while the percentage of geophytes increased along the seainland gradient. C4 photosynthesis was the most prevalent photosynthesis modality in slack areas; most plants had evident Kranz anatomy (Fig. 3a) and were characterised by amphistomatic leaves (Fig. 3b) and abundant aerenchyma in the roots (Fig. 3c). The amount of phenolic compounds present at the root level (Fig. 3d, e) was the lowest in the frontal dunes (Fig. 3f). For continuous traits, the frontal dunes and backdunes were significantly different. The thickness of the adaxial cuticle was higher in the backdunes (Fig. 3g), while limb thickness and LDMC showed the highest and lowest values, respectively, in the frontal dunes (Table 4; Fig.  3h, i). The remaining traits presented a mixed pattern among the three habitat groups or did not show any significant differences.
The PERMANOVA results showed significant differences at the habitat group level for discrete traits and at the transect level for continuous traits (Table 5). The interaction between habitat group and transect was only highly statistically significant for continuous traits.

Relationships Between Species Traits and Environmental Factors
Species traits, environmental factors, and species abundance were highly correlated (RLQ, P < 0.01). The major variation was captured by the first RLQ axis (56%), reflecting the distribution of the plots along the soil gradient (Fig. 4). The organic matter content was highly correlated with the first axis of RLQ (P = 0.01), while soil moisture content was negatively correlated with RLQ2 (P = 0.01), that explained considerably Fig. 2 Non-metric multidimensional scaling diagram based on the similarity (measured by the Bray-Curtis index) among the plots. Different habitat groups are indicated with different symbols and colours. Abbreviations of the sampling localities: T11-T19, plots 1-9 of transect 1; T21-T29, plots 1-9 of transect 2; T31-T39, plots 1-9 of transect 3. All shown plant taxa have a Pearson correlation coefficient >0.7 with the two axes less of the variance (23%) and ordered the plots along a gradient linked particularly to the soil moisture (see Table S3 for a detailed summary of the soil parameters for each habitat group). There were no significant correlations between the other soil variables and the RLQ axes (Fig. S1).
Regarding species traits, invasive taxa, limb thickness, C4 and CAM photosynthesis, and the presence of trichomes, leaf crystals, and aquiferous parenchyma were significantly associated with RLQ1, while the presence of aerenchyma in the roots was significantly associated with RLQ2 ( Table 6). Some of the species traits that were significantly associated with the RLQ axes were correlated with environmental factors. Interestingly, invasive plants showed a positive correlation with the soil moisture content (Fig. 5). Plots with high organic matter content (positive RLQ1 in Fig. 4) were characterised by species with C4 photosynthesis and a low density of leaf trichomes (Fig. 5). These plots (especially T18 and T38) were part of the backdune and slack (T19) plant communities. Finally, plants that colonised the most humid part of the seainland gradient (temporarily wet dune slacks) were characterised by the presence of root aerenchyma (Fig. 4).

Discussion
In the present study, we highlighted how edaphic variables may affect plant community composition through the selection of key functional traits. Among the 23 functional traits, invasive status, C4 photosynthesis, the presence of trichomes on the leaves, and root aerenchyma showed a significant correlation with several soil parameters, such as soil moisture, organic matter content, and pH (Fig. 5). Soil parameters have been shown to be related to vegetation patterns in other Mediterranean coastal dunes and salt marshes (Fenu et al. 2012(Fenu et al. , 2013Ruocco et al. 2014;Landi and Angiolini 2015).
Plants from the frontal dunes were characterised by greater limb thickness as well as lower adaxial cuticle thickness and LDMC than those found in the backdunes. These results could be linked to the presence of aquifer parenchyma in the leaves, which we found in Cakile maritima, Salsola kali, Convolvulus soldanella, Echinophora spinosa, and Eryngium maritimum, which are all target species of Italian beaches and dunes (Angiolini et al. 2018). It seems that the characteristic of succulence is strongly associated with a more saline environment (Mahdavi and Bergmeier 2016), in which case a thin cuticle provides sufficient protection. In contrast, only the Poaceae -Calamagrostis arenaria, Sporobolus pumilus, and Thinopyrum junceum-in the frontal dunes showed sclerophyllic characteristics and possessed a thick cuticle and highly protected stomata to limit water loss (Dickison 2000;Koch et al. 2009) (Fig. S2). As seen in Ciccarelli (2015), the LDMC values of psammophytes were lower in the upper beach and in embryonic dunes, and they increased along the shoreline-inland gradient where the environmental conditions became less harsh, confirming our hypothesis that ruderal species are dominant. The dominant species in the backdunes were either annuals, such as Festuca fasciculata, which escapes the period of greatest water stress (Dickison 2000), or hemicryptophytes, such as Lomelosia rutifolia. These plants show different   Table 1 characterised by plants that differentiated resistant tissues such as sclerenchyma in the leaves and a woody and resistant underground system, as seen in Fumana procumbens, or a parenchymatous underground system as in Lomelosia rutifolia, Silene otites (with calcium oxalate crystals), and Seseli tortuosum (with starch) (Fig. S2). These results partly confirm our hypothesis of the dominance of stress-tolerant species in backdune communities because this habitat group was a heterogeneous environment formed by both herbaceous communities of annuals and shrub communities of perennials showing a mix of ruderal and stress-tolerant traits. Slack areas were characterised by high values of organic matter and soil moisture contents (Table S3) and were dominated by plants that adopted C4 photosynthesis with amphistomatic leaves and abundant aerenchyma in the roots (e.g. Schoenus nigricans, Scirpoides holoschoenus, and Tripidium ravennae). The dominant species in the slack areas, Scirpoides holoschoenus, showed xeromorphic traits (abundant sclerenchyma, thick cuticles, and high LDMC values), as well as typical hydrophytic traits (aerenchyma in the roots and unprotected stomata). There is some functional similarity between backdunes and slack areas that could be explained by the fact that slack areas are generally wet for only a few months; in fact, they experience the same arid conditions as backdunes for most of the year, confirming our hypothesis that plants of slack areas show adaptations to both flooded and arid environments. Similar results have also been found in coastal plants of wet dune slack regions along the Atlantic coast (Bona et al. 2020). The presence of stomata on both the adaxial and abaxial surfaces of the leaves could be interpreted as an adaptation to sand burial, which is a frequent event in maritime sand dune ecosystems, and the verticalisation of the leaves is associated with high light and low water availability (Melo-Júnior and Boeger 2017). Moreover, plants with amphistomatic leaves can increase their photosynthetic efficiency to develop more rapidly (Muir 2015) and to better contend against sand burial (Ciccarelli et al. 2009). In addition, in plants that experience periodic flooding (slack areas), the amphistomatic leaf can act as an advantage for plants living in full-sun environments because it increases their photosynthetic capacity at a time of rapid water fluctuation in the soil (Mott et al. 1982;Koch et al. 2009). In slack areas, most of the perennial herbaceous plants (67%) had aerenchyma in their roots. Also, this feature could be linked to the periodic rainwater accumulation in the interdunal regions, which Exoderm thickness 0.08 −0.01 Fig. 5 Results of fourth-corner analysis of the relationships between plant functional traits and soil variables. Grey boxes represent non-significant correlations (P > 0.05); red boxes represent positive and significant correlations (P < 0.05); blue boxes represent negative and significant correlations (P < 0.05). Functional trait abbreviations are given in Table 1 causes frequent decreases in soil oxygen. Aerenchyma is a tissue with large intercellular spaces that facilitate the circulation of oxygen in the roots to maintain aerobic respiration (Brändle and Crawford 1987;Grosse and Mevi-Schulz 1987;Dickison 2000). Moreover, it seems that increased oxygen supply by aerenchyma outweighs resistance against the mechanical strength of roots in salt marsh vegetation (Veldhius et al. 2019).
Regarding life form, we observed an increase in geophytes along the shoreline-inland gradient. This finding is not entirely consistent with the data of Torca et al. (2019), who analysed Atlantic European coastal dunes. The difference may be linked to the different climatic conditions: Mediterranean coastal dunes are characterised by higher levels of aridity than Atlantic coastal systems. In fact, these plants showed resistance structures such as rhizomes, bulbs, and storage roots and stems that hold water and starch reserves to protect the bud and promote regrowth during favourable periods when rainfall is higher. Another reason that could explain the increase in geophytes along the shoreline-inland gradient is the presence of slack areas in the inner part of the dunes which host many species with this life forms. Slack areas have not been investigated by Torca et al. (2019), who only studied mobile and fixed dunes.
Finally, the plots characterised by high levels of soil moisture were more likely to be colonised by invasive plants. The present results are in accordance with the recent findings of Tordoni et al. (2019), who suggested that native species in coastal sand dunes display functional traits associated with stress resistance, whereas alien species are characterised by traits linked to higher resource acquisition capacity. Nevertheless, it is worth pointing out that the invasion process is strongly context-dependent, and to reach an exhaustive description of the process, it would have been necessary to investigate how soil moisture content varies throughout the year.
In the literature, it has been reported that in embryonic and mobile dunes, functional diversity is lower than that in fixed dunes (Ricotta et al. 2012;Torca et al. 2019). This result may be linked to the high specialisation of plants found in the frontal dunes, which are well adapted to this extreme environment. In fact, environmental filtering is suggested to increase functional similarity (Kembel and Hubbell 2006), while trait divergence is expected under low environmental constraints (Conti et al. 2017). In our study, the PERMANOVA analyses of the CWM values showed that if we consider only discrete traits, the variation can be significantly explained at the habitat group level, while continuous traits were useful to significantly explain the variation at the transect level. Similarly, Mahdavi and Bergmeier (2016) studied coastal sand habitats across different biogeographical regions using discrete traits and found that each habitat type had a similar functional composition across different regions, while different sand habitats in one region had fewer functional groups. Our results highlighted that the different nature of the functional traits used in this study could explain the assembly processes in plant communities and their ecological functions at different spatial scales, from the fine scale of the habitat group level to the coarser scale of the transect level.
Regarding our last goal of using analysis of the relationship between plant functional traits and environmental factors to provide a basis for management strategies, our results highlighted the importance to include root traits in functional studies. While fine roots play an important role in ecosystem functioning, these characters are underrepresented in global trait databases (Iversen et al. 2017). Recently, it has been shown that intraspecific variability in root density can influence sediment stability in salt marshes, highlighting the importance of incorporating easy-to-measure root traits into coastal vegetation monitoring programs (De Battisti et al. 2019).

Conclusions
To the best of our knowledge, this study is one of the first to include measurements of many anatomical traits related to the leaves, stems, and roots of coastal dune ecosystem species. In summary, our results suggest the following: (1) The analysis based on the CWM values is consistent with the ecological constraints along the shoreline-inland gradient improving knowledge of coastal dune plants especially from an anatomical point of view. The frontal dunes were characterised by the presence of ruderal communities who invested in traits (such as succulent leaves, as evidenced by the high limb thickness and low LDMC values) that are well adapted to the harsh environmental conditions of this habitat group. Farther from the sea, the backdune vegetation was characterised by a mosaic of annual herbaceous communities and perennial shrubs showing a mix of ruderal and stress-tolerant traits (clonality, sclerified leaves, high LDMC values, and root phenolics) that are consistent with the less-extreme ecological conditions. Finally, the slack areas seemed to be considerably interesting from a functional point of view, showing adaptations to both flooded and arid environments. They were dominated by plants that had adopted C4 photosynthesis and had both amphistomatic leaves and abundant aerenchyma in the roots.
(2) The RLQ analysis highlighted strong relationships among species abundance, trait occurrence, and soil variables. The invasive status, C4 photosynthesis, leaf trichomes, and presence of aerenchyma in the roots were significantly correlated with soil moisture, organic matter content, and pH. These results seem particularly important for us because they demonstrate the usefulness of anatomical traits (especially those of the root system) in understanding the functional strategy adopted by plants.
(3) Invasive species tended to occupy plots with high levels of soil moisture content and were not abundant in the habitat groups with more arid conditions. Increasing the spatial extent of the study area, investigating how soil moisture content varies throughout the year, and integrating other functional traits, such as ecophysiological or regenerative characteristics, into the study may allow for the development of a more comprehensive functional framework of the invasion process. In addition, we suggest including information regarding root systems (such as dry matter content, aerenchyma presence, root reserves) into coastal monitoring programs because they are directly linked to soil parameters-such as sediment stability or moisture contentuseful in coastal dune management and protection.
Funding Open access funding provided by Università di Pisa within the CRUI-CARE Agreement.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.