Effect of Invasive Alien Species on the Co-Occurrence Patterns of Bryophytes and Vascular Plant Species—The Case of a Mediterranean Disturbed Sandy Coast

Cross-taxon analyses can explain patterns of interaction between taxa and their application in conservation studies can drive management actions. In a coastal sand dune system characterized by a high human pressure, we explored the co-occurrence patterns between vascular plants and bryophytes, with a focus on how the occurrence of invasive alien species (IAS) can affect those taxa and their relationships. Species congruences were evaluated at the community level considering taxonomic and functional diversities. Predictive co-correspondence analysis (Co-CA) was applied to quantify the strength of vascular plant communities in predicting bryophytes species composition. The relationship between the composition of vascular plants and bryophytes was significant, even if weak. Altitude and percentage of bare soil cover are the environmental variables exerting greater influence on the two taxa. The presence of IAS affects communities in an opposite way: for vascular plants, species richness increases with the presence of invasive alien species; for bryophytes, IAS’s presence has a low but significant negative influence, both on species richness and in terms of functional diversity. Results give elements for future studies on the effect of IAS on the bryophytes colonizing coastal sand dunes.


Introduction
Coastal environments represent transitional habitats, characterized by abiotic factors that define a strong gradient of vegetation successions along the sea-land direction [1][2][3]. In terms of biodiversity and landscape composition, their high ecological value is expressed by dynamic interactions between different plant communities in a very restricted area [4,5]. In these ecosystems, an essential ecological function is related to biological soil crusts (BSCs), i.e., ecological pioneer communities of living organisms present on the soil surface in arid and semiarid ecosystems. Due to their ability to fix sand, BSCs promote the increase of organic matter content in the soils, facilitating the existence of the suitable conditions for the establishment, colonization and growth of vascular plants [6,7]. Among biological soil crusts (BSCs), moss crusts represent an advanced stage of development [6]. Regardless of this fundamental ecological part played by mosses, ecological investigations of coastal ecosystems have been mainly focused on vascular plant communities, rarely involving cryptogamic taxa [4]. The few comparative studies on the coastal Mediterranean bryophytic flora revealed that the species pool is poor and characterized by little variability; communities are composed by few species showing acrocarpic growth forms that colonize xeric, sunny and often disturbed habitats, like the arenicolous cosmopolitan Tortella flavovirens. Compared to the Atlantic coastal dunes, Mediterranean ones host few pleurocarpus species, which are more demanding in terms of water and soil availability [4,8].
The lack of knowledge on this crucial, but neglected, taxon could be filled using a cross-taxon approach, i.e., studying the congruence among different biological groups and using the inferred relationship to predict compositional patterns in unknown environments [9]. A significant number of studies on cross-taxon congruence show that terrestrial species richness can be predicted by vascular plant richness and bio-indication [10]. Despite this information, few examples deal with dune ecosystems [11]. Cross-taxon studies exploring covariation patterns between vascular plants and bryophytes are scarce and mainly focused on forest ecosystems [12][13][14][15][16] and nival and riverine habitats [17], showing contrasting results [14].
In the Mediterranean Basin, well preserved coastal habitats became rarer as the anthropogenic pressure increased, mainly consisting in touristic settlements, habitat fragmentation and inappropriate land management, such as reforestation with alien species [11,[18][19][20][21]. Invasive alien species (IAS) are recognized as an important threat to biodiversity worldwide [22] and especially in rare and extreme environments (such as alpine zone or isolated islands, [23][24][25]). On dune ecosystems, in particular, IAS alter the structure of native communities in different ways: reducing functional, phylogenetic and taxonomic diversities and excluding closely related and ecologically similar taxa [26][27][28][29][30].
Despite the well-known urgency to preserve the fragility of coastal ecosystems from IAS in the Mediterranean Basin [2,31,32], the literature on cross-taxon relation involving IAS is still limited. In coastal areas where Acacia spp, Eucaliptus spp and Pinus spp were used for reforestation, those alien species are limiting light availability to bryophytes and photophyles of lichens, but they also constitute a long-degrading litter, affecting negatively soil availability [11].
Moreover, notwithstanding that in recent decades plant functional trait-based approaches in ecological studies grew enormously [33], also focusing on dunes and biological invasions [34][35][36], we noticed that very few studies investigated the relationship among functional diversity of vascular plants and bryophytes [37].
Increasing urban sprawl and littoralization on coastal areas are squeezing dune ecosystems between coastal erosion and the expanding human pressure: in the last 50 years, those habitats have been increasingly disturbed and threatened by several hazards, by depleting the biodiversity and transforming their structure [2,5,19,21]. This situation is posing questions on the appropriate management to preserve the remnants of coastal dune ecosystems and the need to focus on the role of invasive alien species, also considering the recent legislative instruments that Europe provided for IAS control and eradication (EU Regulation 1143/2014 on Invasive Alien Species).
In this context, we investigated a Mediterranean coastal dune characterized by a high human pressure, focusing on vascular plants and bryophytes composition to learn how those disturbed ecosystems respond to environmental features and to the presence of IAS, giving useful hints for their management.
Based on a probabilistic sampling, we estimated the congruence patterns (both in terms of richness and composition) between vascular plants and bryophytes. Furthermore, considering the high impact of IAS in the study area, we are interested in analyzing the possible negative influence of IAS on these relationships. In detail, our purpose is to answer to the following research questions: (i) can the composition of vascular plants predict the composition of the bryophytes communities? (ii) What are the environmental factors related with the covariation of the two communities? Which are the factors better explaining the bryophytes diversity pattern? (iii) Does the presence of IAS have an impact on vascular plants and bryophytes communities and on their functional diversity?

Study Area
The study area is located in Villasimius, Southern Sardinia, Italy. The site is the coastal sand dune of "Porto Giunco" within the Site of Community Interest ITB040020 "Isola dei Cavoli, Serpentara, Punta Molentis", fully included in the marine protected Area Capo Carbonara (Figure 1). Lithologies are essentially referred to the batolite of Sarrabus emerged at the end of the Ercinic orogenesis in the Middle Carboniferous, mainly composed of granites and granodiorites subjected to tectonic movement during the alpine orogenesis [38]. The bioclimate for the area is described as a Lower Thermomediterranean, Lower Dry, Euoceanic Strong, with a mean temperature of 17.8 • C and 561 mm of annual rain [39].

Study Area
The study area is located in Villasimius, Southern Sardinia, Italy. The site is the coastal sand dune of "Porto Giunco" within the Site of Community Interest ITB040020 "Isola dei Cavoli, Serpentara, Punta Molentis", fully included in the marine protected Area Capo Carbonara (Figure 1). Lithologies are essentially referred to the batolite of Sarrabus emerged at the end of the Ercinic orogenesis in the Middle Carboniferous, mainly composed of granites and granodiorites subjected to tectonic movement during the alpine orogenesis [38]. The bioclimate for the area is described as a Lower Thermomediterranean, Lower Dry, Euoceanic Strong, with a mean temperature of 17.8 °C and 561 mm of annual rain [39]. The dominant vegetation can be described as Sardinian psammophilous geosigmetum of coastal sand dunes (Cakiletea, Ammophiletea, Crucianellion maritimae, Malcolmietalia, Juniperion turbinatae) and by the presence of small Juniperus sppl communities composed by J. macrocarpa and J. phoenicia var. turbinata (habitat 2250* Coastal dunes with Juniperus spp., [40]). Dominant species are Juniperus oxycedrus ssp. macrocarpa and Pistacia lentiscus, accompanied by Pancratium maritimum, Lotus cytisoides and Crucianella maritima. Remnants of habitats of high naturalistic value are strongly menaced by anthropic pressure, namely urbanization, touristic exploitation and unadvised reforestation that were carried out using alien species, especially in the backdunes, mostly using Carpobrotus spp., Agave spp., Eucalyptus spp. and Acacia spp. Due to this situation, in recent decades the area was involved in different naturalistic engineering projects in order to protect, conserve and restore the typical dune vegetation. The projects included the positioning of cane structures for sand retention and the consequent recovery of the dune profile. The eradication of invasive species (Carpobrotus spp., Agave spp. and Acacia spp.) in the priority habitats were accompanied with the use of naturalistic engineering techniques, such as the positioning of coconut fiber bio-nets to avoid erosion, particularly subsequently to eradication activities. After the eradications, the sites were sown with seeds and native plants multiplied in nurseries using germplasm previously collected on site [41]. The dominant vegetation can be described as Sardinian psammophilous geosigmetum of coastal sand dunes (Cakiletea, Ammophiletea, Crucianellion maritimae, Malcolmietalia, Juniperion turbinatae) and by the presence of small Juniperus sppl communities composed by J. macrocarpa and J. phoenicia var. turbinata (habitat 2250* Coastal dunes with Juniperus spp., [40]). Dominant species are Juniperus oxycedrus ssp. macrocarpa and Pistacia lentiscus, accompanied by Pancratium maritimum, Lotus cytisoides and Crucianella maritima. Remnants of habitats of high naturalistic value are strongly menaced by anthropic pressure, namely urbanization, touristic exploitation and unadvised reforestation that were carried out using alien species, especially in the backdunes, mostly using Carpobrotus spp., Agave spp., Eucalyptus spp. and Acacia spp. Due to this situation, in recent decades the area was involved in different naturalistic engineering projects in order to protect, conserve and restore the typical dune vegetation. The projects included the positioning of cane structures for sand retention and the consequent recovery of the dune profile. The eradication of invasive species (Carpobrotus spp., Agave spp. and Acacia spp.) in the priority habitats were accompanied with the use of naturalistic engineering techniques, such as the positioning of coconut fiber bio-nets to avoid erosion, particularly subsequently to eradication activities. After the eradications, the sites were sown with seeds and native plants multiplied in nurseries using germplasm previously collected on site [41].

Sampling Design and Data Collection
We delimited the main polygons of residual dune vegetation located along the 3 km long sand dune of Porto Giunco ( Figure 1) and, to avoid the edge effect, we applied a 10 m internal buffer to the polygons. On the remaining area of about 10 ha (7 polygons), according to a simple random sampling, we selected one hundred 60 × 60 cm squared plots, a standard dimension to sample Mediterranean bryophytes [42][43][44]. In each plot we measured the presence and the cover abundance (in percentage) of vascular plants and bryophytes.
As environmental predictors, we used the plot spatial coordinates, gathered in the field using a GPS (Garmin GPSMAP 64S, UTM WGS84), along with the plot distance from the sea line and the altitude a.s.l. Moreover, each plot was assigned to one of the main four identified vegetation types (bare soil, low scattered Mediterranean scrub, low Mediterranean scrub and tall Mediterranean scrub, Table 1), according to the dominant physiognomy (covering more than 50% of the plot). Vascular plants nomenclature follows Pignatti [45]; vascular plant species were assigned the IAS status following Galasso et al. [46]. Nomenclature for bryophytes followed Hodgetts [47].

Data Analysis
Cross-taxon congruence between vascular plants and bryophytes, their relationships with the environmental predictors and the effect of the presence of invasive vascular plants were assessed using different statistical methods. Firstly, two plot-to-plot compositional distance matrices were calculated separately for vascular plants and bryophytes using the Bray-Curtis dissimilarity coefficient on transformed (log+1) abundance data. A third Bray-Curtis dissimilarity matrix was also calculated, for the same plots, but it considered solely the IAS composition for vascular plant composition, adding a dummy species (0.001 abundance) to all samples to avoid calculation problems where no IAS have been detected. Conversely, no bryophytes IAS were sampled in any plot. The same correlation analyses were also performed downscaling the species abundances to presence/absence and using the Jaccard dissimilarity coefficient.
Species functional traits for both plant communities and bryophytes were used to calculate pairwise functional distances among plots using the package "adiv" [57]. Specifically, considering the different nature of functional traits available (quantitative, ordinal, nominal), we computed two functional dissimilarity matrices (one for vascular plants and one for bryophytes) using a modified version of the Gower distance allowing for the treatment of various statistical types of variables (for more details about the methodology see [58]). To investigate the strength of the correlations among compositional and functional diversity of vascular plants vs. bryophytes and the correlations among compositional and functional diversity of invasive alien vascular plant species vs. bryophytes, we used the r Pearson's correlation coefficients testing their significance with Mantel tests (based on the null distribution of the r Pearson' correlation coefficient on 9999 permutations, selected α level = 0.05). Partial Mantel tests were also performed testing the same relationships described above but using a conditional distance matrix based on the Euclidean distances of the plot variables (i.e., coordinates, altitude, distance from the sea line, vegetation physiognomy).
Furthermore, predictive co-correspondence analysis (Co-CA, see [59,60] for a full description of this method), was used to quantify the strength of vascular plant community data in predicting bryophytes species composition. A second predictive Co-CA was also performed to explore if and how IAS vascular plant composition influences (positively or negatively) the composition of bryophytes. A leave-one-out cross-validatory fit percentage was estimated, for each Co-CA, to select the minimal adequate predictive ordination axes. Cross-validatory fit > 0 implicitly validates the model, indicating that prediction is better than the one obtained under the null model [60]. Hence, a permutation test for predictive Co-CA models were performed to assess the significance of each ordination axis (statistical significance level selected α = 0.05, 999 permutations); in addition, environmental variables were fitted onto Co-CA ordination axes aiming at understanding the meaning of axes in term of the underlying ecological gradients. The "envfit.coca" function in the "vegan" R package [61] was used at this scope and the associated goodness of statistic fit was the squared Pearson correlation coefficient (r 2 ).
Finally, using generalized linear models (GLM) we tested for linear relationships of biotic (IAS presence) and abiotic predictors on (1) species richness and (2) functional diversity of bryophytes in the study area. The selected measure of bryophytes functional diversity for the GLM was the plot-based Rao's Q proposed by Ricotta et al. [62] (all functional traits available for bryophytes were used for this calculation). As predictors, in both models, variables listed in Table 1 were used. A step-wise variable selection procedure was adopted to consolidate the Minimal Adequate Model (MAM). Akaike Information Criterion (AIC) was selected to find, among all the subsets of reduced models, the final MAM.
Mantel and partial Mantel tests were computed using R package "vegan" [61]; Co-CA and related diagnostics with R package "cocorresp" [63]. All the other statistical analyses were performed using basic R functions like "glm" [64]. For more details on the GLM specifications and selected options, please, see Supplementary Materials I.

Results
Eleven bryophytes and 55 vascular plants were detected. Five vascular plants were classified as invasive alien species (IAS): Acacia saligna, Agave americana, Carpobrotus spp., Myoporus tetrandum and Opuntia ficus-indica (Table 2). IAS occurred in 13% of the sampled plots: Carpobrotus spp. was the most abundant one and it was recorded in 8% of the sampled plots with a variable coverage, ranging from 25% to 80%.
Among Bryophytes, Tortella flavovirens was the most common species. This species, that could be classically described as a pioneer species, was the only bryophyte present in plots where IAS were recorded. The mean IAS coverage per plot was less than 5% (Table 1). The correlation between vascular plants and bryophytes species richness was positive but not significant (r = 0.15; p > 0.05). However, predictive Co-CA on the species composition for the two groups showed that vascular plant communities have a certain predictive power on the composition of bryophytes (plants explain 21.4% of the cumulated variance in bryophytes when significant axes are considered, Figure 2 and Table 3). In the Co-CA ordination space, altitude and percentage of bare soil were significantly related with the compositional differentiation of both vascular plants and bryophytes communities (r 2 = 0.15, p < 0.05 for altitude, r 2 = 0.149, p < 0.01 for bare soil, Figure 2). This result was also corroborated by GLMs (Tables 4 and 5, Supplementary Materials I): percentage of bare soil and vegetation type were the main variables statistically supported for bryophytic species richness (Table 4) and functional diversity (Table 5).   Decoupling the two communities, we found a low but significant correlation among the two taxa, using both presence/absence data (Jaccard dissimilarity, Mantel test r = 0.2, p < 0.001) or abundance data (Bray-Curtis dissimilarity, r = 0.15, p < 0.001, Figure 3).   The correlation remains significant also when we control for environmental variables via Partial Mantel tests (geographic coordinates, distance from the sea, altitude and physiognomy of the vegetation. Jaccard dissimilarity r = 0.2 p < 0.001, Bray-Curtis dissimilarity r = 0.13 p < 0.001).
When we investigated the effect of IAS on native communities, we observed that these influenced both the composition and functional diversity of native vascular plants and bryophytes (Figure 4). The predictive Co-CA analysis on bryophytes based on IAS explains only~6% of their variance ( Figure 4, Table 3).   As in the previous analysis, altitude and percentage of bare soil were significantly correlated with the compositional differentiation of both IAS and bryophytes communities (r 2 = 0.119, p < 0.01 for altitude, r 2 = 0.188, p < 0.001, for bare soil, Figure 4).
While vascular species show an increase in the presence of alien species, bryophytes are negatively affected by IAS composition and functional diversity ( Figure 4). IAS presence/absence affects negatively bryophytic species richness (Table 4), while functional diversity of bryophytes is negatively influenced by IAS presence (Table 5).

Discussion
In the Mediterranean basin, coastal ecosystems have been widely investigated to understand patterns of species distribution along natural gradients [4,65,66]. In accordance with previous studies [4], we found that vegetation is distributed along a sea-inland gradient (i.e., altitude and bare soil): both vascular and bryophytes communities follow this gradient, pointing out the pioneer role of bryophytes in establishing soil, trapping and holding moisture [7]. In particular, bryophytes associated to abundant bare soil cover are Barbula unguiculata (Pottiales) and Bryum dichotomum (Bryales), stress tolerant species living closer to the sea, typical of the communities described as ground pioneer communities with postpioneers, acidiclines with neutroclines, with xerophilic tendency (Barbuletea unguiculatae Mohan 1978) able to colonize areas where sea aerosol and winds are stronger [8,53]; another indicative species is Ptychostomum imbricatulum, known to be an indicator of moderate/high anthropic impact (Meso-Euhem) [8,53].
In the same way, among vascular plants we can relate the abundance of bare soil with the colonizing gradient described by the presence of Salicornia ramosissima, a salt tolerant plant, typical of annual communities, Crucianella maritima and Pancratium maritimum typical of Western Mediterranean, chamaephytic, hind-dune communities that develop on semipermanent dunes (Habitats Directive, 2210 Crucianellion maritimae fixed beach dunes) and Eryngium maritimum, a species that grows typically on sand and shingle beaches, foredunes and yellow dunes, as well as in semifixed grey dunes. Based on our results we can affirm that a weak relation between bryophytes and vascular plant exists along this gradient, hence the vascular plant community on Mediterranean sandy dunes can be used as a predictor of bryophyte community. This result can be interpreted as linked to the response to the abiotic factors affecting species richness and functional diversity of bryophytes in the study area: percentage of bare soil and vegetation type were the main variables statistically supported for bryophytic species richness and functional diversity.
The most common bryophyte is Tortella flavovirens: This pioneer species is frequent in Mediterranean coasts and it usually forms stable cushions leading to soil consolidation (see also [4]) and it is usually associated with foredunes and retrodunal environments. In our study it is the only bryophyte that co-occurs with IAS.
Regarding IAS presence and effect on the vascular plants and bryophytes communities, we observed an opposite response: the correlation for bryophytes is negative in terms of composition and functional traits, while it is positive for vascular plant composition. This relationship is not linked to the response to environmental factors and gives us an important element for future studies on the effect of IAS management on the bryophytes colonizing the coastal sand dunes.
All IAS recorded in our investigation are perennial species: vascular plants might benefit from the IAS presence because they might increment soil and water availability all year long. The negative effect on bryophytes species due to alien species could be associated to the fact that, on dunes, many bryophytes grow under shrubs and trees, taking advantage of shade and a greater water availability [38]. In Mediterranean natural communities, this role is usually covered by bushes or trees such as Juniperus spp., which grow at a slower rate if compared with IAS [4], suggesting that invasive alien species substituting native species can induce an indirect reduction on the diversity of bryophytes. Within IAS, Carpobrotus spp. grow closer to the sea and its catastrophic effect in reducing functional and taxonomic diversity of vascular plants is well investigated [67,68]. This fast-growing horticultural succulent species from South Africa affects negatively the presence of native species limiting space availability and its roots cause soil acidification, leaving a long-lasting effect even after its removal.
Nevertheless, co-occurrence among IAS vs. bryophytes and the eventual disruption of the relationship observed among bryophytes and vascular plants was not detectable since IAS and bryophytes coexist only in the 15% of the invaded plots (i.e., Tortella flavovirens with Carpobrotus and Agave americana), definitely data too fragile to extract generalities. We must remember that the area involved in our study is highly compromised by human impact and the vegetation is not shaped exclusively by natural factors. Our results must consider that the high anthropic pressure on natural communities could be reinforced by the presence of IAS, which were mostly intentionally introduced or escaped from private gardens [69,70]. IAS that were traditionally planted to stabilize sand dunes such as Acacia saligna, might induce a contrary effect in the long period, affecting and disrupting not only native vascular plants communities [71] but also the bryophyte ones. However, to investigate how the presence of IAS influences ecological successions through time, long-term and specific design (e.g., BARCI) and analyses will be necessary.
In spite of the relevance in explaining ecological invasion by IAS, a cross-taxon approach is still rare and results are still contrasting. We found that vascular plants and bryophyte communities are both influenced by altitude and bare soil. However, the two communities respond differently to the presence of IAS. Indeed, vascular plant community richness is promoted by IAS presence perhaps because the latter promote greater availability of water or because it follows the rich get richer hypothesis [72], contrasting the invasion paradox [73] and supporting recent global meta-analyses [74]. On the contrary, the presence of IAS has a negative and low but significant influence on bryophytes in terms of species richness and functional diversity. In this study, the low co-occurrence IAS vs. bryophytes does not allow us to gain strong conclusions on the IAS's effects on cross-taxon relation. Future studies should focus on the influence of spatial scale, introducing a multiscale approach, to deepen the functional approach and investigate the pedological component to understand invasion patterns in very fragile ecosystems such as Mediterranean coasts.