Pollination networks along the sea-inland gradient reveal landscape patterns of keystone plant species

Linking the functional role of plants and pollinators in pollination networks to ecosystem functioning and resistance to perturbations can represent a valuable knowledge to implement sound conservation and monitoring programs. The aim of this study was to assess the resistance of pollination networks in coastal dune systems and to test whether pollination interactions have an explicit spatial configuration and whether this affect network resistance. To this aim, we placed six permanent 10 m-wide belt transects. Within each transect we placed five plots of 2 m x 2 m, in order to catch the different plant communities along the dune sequence. We monitored pollination interactions between plants and pollinators every 15 days during the overall flowering season. The resulting networks of pollination interactions showed a relatively low degree of resistance. However, they had a clear spatial configuration, with plant species differently contributing to the resistance of pollination networks occurring non-randomly from the seashore inland. Our results evidenced that beside contributing to the creation and maintenance of dune ridges, thereby protecting inland communities from environmental disturbance, plant species of drift line and shifting dune communities have also a crucial function in conferring resistance to coastal dune pollination networks.

Scientific REPORTS | (2018) 8:15221 | DOI: 10.1038/s41598-018-33652-z interconnected through flows of energy and materials at the landscape scale, and many examples of plant communities which naturally form continuous patterns in the landscape in response to environmental gradients have been described (e.g., river and lake edges, salt marshes).
Coastal dune systems are among the most outstanding examples worldwide of spatially connected plant communities 17,18 . Environmental disturbance in the form of salt spray, wind, soil movement and wave regime shapes coastal dunes selecting species according to their tolerance to the environmental disturbance and creating a precise sequence of ecologically distinct plant communities from the seashore inland, the so-called coastal zonation 19,20 . The coastal zonation reflects the progressively more sheltered environments that are created with increasing distance from the shoreline. Pioneer salt spray and sand blasting tolerant plant communities developing nearby the seashore are responsible for the formation and the building-up of dune ridges 21 that act as a physical barrier to different stressors (e.g., wind, sand-blasting, waves), ultimately creating suitable conditions for the establishment of more disturbance-sensitive communities typical of semi-fixed and fixed dunes 22 .
Although the major dune-forming species are wind-pollinated grasses, animal-pollinated species occur in all plant communities of coastal dune systems and concur to the development of dune ridges (e.g., Cakile maritima Scop., Echinophora spinosa L., etc. 19 ). Many of them are listed as diagnostic species for habitats definition in the Interpretation Manual of European Union Habitats (EU 28) 23 and they are used as indicator species to define habitats' conservation status 17,24 . Moreover, nowadays they are facing a huge decline because of coastal dune degradation and surface contraction 21,25 , and many of them have become a topic of increasing conservation concern 26,27 .
In the light of these considerations we assessed the function of pollination in coastal dune systems by answering the following questions: (i) does the richness of animal-pollinated plant and pollinator species show non-random spatial patterns from the seashore inland? (ii) Are coastal dune pollination networks resistant to perturbations (i.e., nested, modular and redundant) during the overall flowering season? (iii) Do plant species growing at different distances from the sea differently contribute to the resistance of coastal dune pollination networks?

Material and Methods
Study area. We selected an area along the North Adriatic coast (Italy), in the peninsula of Cavallino (45°26′22.14″N; 12°26′58.71″E), where a complete sequence of coastal dune communities, from the drift line to the inland fixed dune, can be observed. Sites consist of narrow, recent dunes (Holocene), that generally occupy a narrow strip along the seashore, bordered by river mouths and tidal inlets, mostly fixed by docks. The study area has a Temperate Oceanic bioclimate 28 . The annual average temperature is about 13 °C and the annual average precipitation is 831.5 mm 22 . Sediments are predominantly sandy carbonate deposits of dolomitic origin transported by river Piave and Tagliamento. The study area is under protection and has been subjected to conservation plans aimed at reducing natural erosion and human trampling 29 .
In the study area, plant community sequence starts with the pioneer, nitrophilous community dominated by annuals of the drift line zone (C. maritima plant community). Being exposed to wave inundation, salt spray and wind stress, the community is often patchy and fragmented. Landwards, the C. maritima plant community is followed by the shifting dune community dominated by dune-forming plants such as Elytrigia juncea (L.) Nevski and Ammophila arenaria subsp. australis (Mabille) Laínz. A. arenaria subsp. australis is the dominant species and is responsible for stabilizing and building up the foredune by capturing blown sand and binding it together with its tough, fibrous rhizome system 30 . Inland of the foredune crest, protection from physical stresses increases and the vegetation evolves towards more structured forms. The coastal sequence includes perennial xerophilous grasslands of the semi-fixed dunes dominated by dwarf shrubs (e.g., Fumana procumbens (Dunal) Gren. & Godr., Thymus pulegioides L.), lichens (e.g., Cladonia foliacea (Huds.) Willd., C. furcate (Huds.) Schrad.) and mosses (e.g., Syntrichia ruraliformis (Besch.) Cardot), and interdunal depressions where a community of Erianthus ravennae (L.) H. Scholz and Schoenus nigricans L. can be observed. Lastly, the sequence ends with the xerophilous shrublands with Erica carnea L. and Osyris alba L. and the xerophilous woodlands with Quercus ilex L., Pinus pinea L. and P. pinaster Aiton of fixed dunes 31,32 (Fig. 1).

Data collection.
At the beginning of the season we placed 6 permanent belt transects. Transects were 10 m wide and 606 ± 25 m (mean ± SD) long. The distance between transects was set at 50 m. Along each transect we identified plant communities and their spatial extent based on the habitat map of the Veneto region (available at https://www.regione.veneto.it/web/agricoltura-e-foreste/download#IT3250003; scale 1:10.000). Sampled communities were: (i) the pioneer community of the drift line, (ii) the pioneer community of shifting dunes, (iii) the xerophilous grasslands of semi-fixed dunes, (iv) the xerophilous shrublands of fixed dunes and (v) the xerophilous woodlands of fixed dunes (Fig. 1). All plant communities of the coastal system were sampled with the exception of interdunal depressions because of the almost total absence of animal-pollinated species. Within each transect we placed one plot of 2 m × 2 m per plant community, for a total of 5 plots per transect, and 6 replicates per plant community.
Due to the particular characteristics of coastal dune communities (e.g., scantily covering vegetation, prevalence of wind-pollinated species), within each 10 m wide belt transect and the spatial extent of each community, we placed the plots according to a preferential sampling design so as to catch at least one animal-pollinated species. Plots distance from the sea was measured through a trekking GPS (Garmin GPSmap 62 s). In each plot, plant species presence was recorded, and the flowering phenology of animal-pollinated species was monitored every fifteen days, from the beginning (April 03-04-2017) to the end of the flowering season (August 01-08-2017), for a total of nine surveys. The nine surveys were carried out in warm and sunny days 33,34 . Flowering monitoring started at the opening of the first flower and ended when individual plants no longer possessed any flower with anthers 35 . Furthermore, we quantified the length of flowering of each animal-pollinated species by counting the number of surveys in which a given species was observed flowering. During each survey, we recorded visiting pollinators by counting the number of visits to each plant species. Each plot was monitored for 14 minutes, split into two 7-min sets distributed during two daily intervals (from 10 a.m. to 1 p.m. and from 1 p.m. to 4 p.m.) to ensure the observation of pollinators showing different daily periods of activity. We set the starting time of monitoring at 10 a.m. and the finishing time at 4 p.m., because, according to preliminary observations, pollination events were infrequent before 10 a.m. and after 4 p.m. Since literature reported nocturnal moths as potential pollinators of Silene vulgaris (Moench) Garcke 36 , we also conducted nocturnal surveys in order to observe pollination events mediated by nocturnal moths. However, no pollination contacts could be observed. Our methodology is in line with e.g., Lázaro et al. 33 , Hegland and Totland 37 and Fantinato et al. 13 . In fact, with a few exceptions 38 , this time span has been proven to overlap the daily pattern of activity of the majority of pollinator species in temperate and Mediterranean ecosystems. In total we conducted 252 (i.e. 1,764 min) observation periods (126 in the morning and 126 in the afternoon) during the flowering season. We counted as pollinators only visitors landing on the flower, visiting it for more than 1 s, and getting in direct contact with the floral reproductive organs 37  Sampling completeness. Since we sampled plant communities subjected to different environmental conditions, which might possibly lead to differences in pollinator activity, for each plant community we evaluated the sampling completeness of animal-pollinated plant and pollinator species, and of pollination interactions. We estimated the richness of (i) animal-pollinated plants, (ii) pollinators and (iii) pollination interactions with increased sampling effort, pooling data recorded in each plot during the overall flowering season. We computed the accumulation curves using the cumulative number of plots sampled per plant community as the unit of sampling effort and the Chao 2 estimator of asymptotic species richness 39 . We chose Chao 2 estimator, because it is one of the least biased estimates for a small number of samples [40][41][42] . Sampling completeness was calculated by quantifying the percentage of the asymptotic richness detected by the observed one 42 . The accumulation curves and the Chao 2 estimator were calculated with the R-based package Vegan 43 .
Network parameters. For each survey, observed contacts were organized in an abundance matrix in which rows were plant species, columns were insect species, and entries represented the sum of all interactions observed between each plant and pollinator species. We chose to organize pollination interactions in one matrix for survey (for a total of 9 matrices) to avoid the formation of impossible interactions through pollinator sharing between plant species flowering during different surveys (i.e., forbidden links 44 ). For each pollination matrix, we calculated the most used quantitative indices to describe the structure and resistance of weighted ecological interaction networks at the level of the entire dune system as well as at single species level. At the dune system level we calculated network selectiveness (H′ 2 45 ), weighted nestedness (wNODF 46 ) and quantitative modularity (Q 47 ). A low degree of network selectiveness is linked to high interaction redundancy, which confers resistance to pollination networks by acting as a buffer against species loss. A nested pattern of interactions confers resistance to pollination networks by preventing the occurrence of secondary extinctions 48 , while a modular organization of interactions prevents the propagation of secondary extinctions through the pollination network 49 . The significance of observed results was tested by constructing 1000 randomized networks with identical margin totals as the empirical networks and comparing the observed and random values using the null model 'r2d' 50 .
At the species level we chose indices revealing the contribution of each plant species to network selectiveness, nestedness and modularity. Especially, we calculated species selectiveness (d′ 45 ), which measures the exclusiveness of a species partner spectrum compared with other species in the network; nested contribution (n i 51 ), which Data analysis. Since the extent of plant communities increased exponentially from the seashore inland, spanning from few meters for the pioneer community of the drift line to hundred meters for the coastal wood community, we log-transformed plot distance from the sea before data analysis.
Spatial variations in the richness of animal-pollinated plant and pollinator species were assessed through linear mixed models (LMMs; R-based package nlme 54 ). We regressed the richness of animal-pollinated species in bloom with respect to a quadratic trend, while the richness of pollinators with respect to a linear trend. Especially, we included plot identity as random factor, the log-transformed distance from the sea as independent variable and the richness of (i) animal-pollinated species in bloom and (ii) pollinators as dependent. We used species richness of animal-pollinated species in bloom and pollinators recorded in each survey as replicates.
Subsequently, we calculated the mean log-transformed distance from the sea of each animal-pollinated species by averaging the log-transformed distance from the sea of the plots in which a given plant species was recorded. We assessed spatial variation of plant species selectiveness (d′), nested contribution (n i ), weighted closeness centrality (wCC), standardized connection (c) and participation values (z) by means of linear mixed models (LMMs). We then regressed each index with respect to a linear trend, with the identity of species as random factor, the mean distance from the sea as independent variable and each single index as dependent. If a given plant species was observed to flower for more than one survey, we used values of indices calculated in each survey as replicates.
Furthermore, the spatial variation of the flowering length of animal-pollinated species was tested through a generalized linear model (GLM) using a Poisson distribution and log as link function with respect to a linear trend, the mean log-transformed distance from the sea as independent variable and the flowering length as dependent.

Results
Overall, we recorded 521 contacts between 18 plants and 59 pollinators. The 18 animal-pollinated species belonged to 13 different families. The majority of plant families were represented by a single species (nine families). The most specious families were the Caryophyllaceae (three species), followed by the Fabaceae, Lamiaceae and Rosaceae (both with two species). Plant species were highly specialized for a single plant community; only few species were found in more than one plant community (e.g., E. spinosa, Calystegia soldanella (L.) Roem. & Schult.), and always with few individuals. The 59 pollinators recorded were all insects and belonged to six orders. Among the pollinator orders the most specious were the Hymenoptera (37 species), followed by the Coleoptera (eight species), Diptera and Lepidoptera (six species, both). Overall, 33 pollinator species were found in only one plant community, 15 species in two plant communities, six species in three plant communities, three species in four plant communities, while only Apis mellifera (Linnaeus, 1758) and Bombus pascuorum (Scopoli, 1763) were found in all plant communities.
Sampling completeness revealed that we detected more than 82.75% of animal-pollinated plants and more than 75.50% of pollinator species in all plant communities (Table 1; Fig. 2). For the richness of pollination interactions (Table 1; Fig. 2), detection varied among plant communities, diminishing from the pioneer community of The richness of animal-pollinated species in bloom and of pollinator species notably changed along the sea-inland gradient. Linear mixed models (LLMs) revealed that the richness of animal-pollinated species in bloom followed a unimodal trend (t-value = −2.738; P = 0.011; Fig. 3a), with the peak of species richness coinciding with the xerophilous grasslands of semi-fixed dunes. On the contrary, pollinator species richness decreased significantly from the seashore inland, peaking in correspondence with the pioneer community of the drift lines (t-value = −3.227; P = 0.003; Fig. 3b).
Observed values of network selectiveness, nestedness and modularity were significantly different from random values throughout the flowering season (Table 2). Nevertheless, values of network selectiveness were always high (always >0.5), and coupled with low values of nestedness (<15; with the exception of the 9 th survey), and low values of modularity (almost always <0.6).
Linear mixed models (LMMs) revealed that species selectiveness (d′) significantly increased from the seashore inland, while weighted closeness centrality (wCC) and nestedness contribution (n i ) significantly decreased ( Table 3   The flowering length of animal-pollinated species decreased significantly from the seashore inland (z-value = −2.249; P = 0.024), with C. maritima being the species flowering for the longest period of time in the studied system (6 surveys, approximately 2.5 months).

Discussion
Coastal dune systems are shaped by strong environmental gradients which determine the co-occurrence, in a relatively small area, of different plant communities. The steep gradient allows the presence of communities characterized by a high ecological diversity in terms of environmental properties as well as of species composition, with different pools of plant and pollinator species showing high specialization for a single plant community. In terms of richness of animal-pollinated species, our results confirm xerophilous grasslands of semi-fixed dunes as the most outstanding community of costal dune systems. This pattern was already evidenced by previous studies 19 which observed an increase in plant species richness following the sea-inland gradient, with a slight decrease in the woody community of fixed dunes. Accordingly, richness of animal-pollinated species followed a unimodal trend, reaching the highest values in correspondence with the xerophilous grasslands of semi-fixed dunes. Being less exposed to the environmental disturbance, xerophilous grasslands normally present a higher plant species richness compared to both the pioneer communities of drift line and shifting dunes 19,20 . Interestingly, pollinator species richness did not follow a comparable trend; pollinator richness significantly decreased from the seashore inland, with a peak in the pioneer community of drift lines, which, however, hosts a low number of animal-pollinated species.
In coastal dune systems, few animal-pollinated species (i.e., C. maritima, C. soldanella, E. spinosa, M. marina) are responsible for the maintenance of pollinator species richness. The observed pattern is in contrast with those observed in other plant communities 55,56 . Numerous empirical studies report a positive relationship between the richness of plant species and that of pollinator species 57,58 . The positive impact of the plant species richness on the richness of pollinator species has been explained by increasing floral resource heterogeneity (nectar and pollen), which increases attractiveness for many pollinator species seeking single and multiple resources 8 . However, it is possible that species with large flowers/inflorescences like C. maritima, C. soldanella, E. spinosa and M. marina, attract more pollinators than plant species with isolated, often poorly attractive flowers/inflorescences like, e.g. Cerastium semidecandrum L., F. procumbens and S. vulgaris, which are typical of the xerophilous communities of the semi-fixed dunes (but see Ebeling et al. 56 ). In this context, the quantification of floral resources (nectar and pollen) produced by each species may provide new insights into the observed mismatch between the richness of plant and pollinator species, thus revealing how few animal-pollinated species can sustain the richness of pollinator species in coastal dune systems.
Another explanation might be related to an unbalanced sampling effort among plant communities. However, our analysis of sampling completeness revealed that the sampling effort allowed us to catch a high proportion of the richness of animal-pollinated plant and pollinator species in all plant communities. A decreasing completeness from the pioneer community of the drift line to the xerophilous woodlands of fixed dunes was detected only for pollination interactions. One possible explanation is that, moving inland, pollination interactions become more selective, resulting in a higher probability to observe single interactions. Chacoff et al. 42 pointed out that, since the Chao 2 estimator quantifies completeness on the basis of infrequent species (or interactions) 40,41 , a difference between the detected number of plant and pollinator species, and that of pollination interactions, could be attributed to the presence of selective interactions, rather than to true undersampling of interactions. However, Survey 1st 2nd 3rd 4th 5th 6th 7th 8th 9th H′ 2 1.000 *** 0.600 *** 0.780 *** 0.592 *** 0.601 *** 0.748 *** 0.801 *** 0.687 *** 0.681 *** WNDOF 0.000 *** 14.102 *** 8.781 *** 9.190 *** 9.018 *** 4.991 *** 11.160 *** 9.484 *** 28.368 *** Q 0.500 *** 0.347 * 0.581 *** 0.521 *** 0.491 *** 0.662 *** 0.162 *** 0.502 *** 0.185 *** Table 2. Network parameters for each survey. H′ 2 : network selectiveness; WNODF: weighted nestedness; Q: quantitative modularity. The significance of observed results was tested by constructing 1000 randomized networks with identical margin totals as the empirical networks and comparing the observed and random values using the null model 'r2d' ( * P < 0.05; *** P < 0.001).  Table 3. Results of the linear mixed models (LMMs) comparing plant species level network indices and mean log-transformed distance from the sea. Each index was regressed with respect to a linear trend, with the identity of species as random factor, the mean distance from the sea as independent variable and each single index as dependent. If a given plant species was observed to flower for more than one survey, we used values of indices calculated in each survey as replicates.
this remains a point that deserves further studies and we cannot exclude that a further sampling effort might have increased the completeness of pollination interactions in the xerophilous communities of semi-fixed and fixed dunes.
Interestingly, in accordance with the decreasing see-inland trend of completeness of pollination interactions, we found that the degree of selectiveness (d′) of plant species increased throughout the coastal dune system, reaching the highest values in the xerophilous woodlands of fixed dunes. High selectiveness entails a high dependence of plant species to pollinator species, thereby resulting in less resistance to perturbations. Conversely, a low plant selectiveness is associated with high pollinator sharing and thus redundancy of interactions. Redundant interactions confer high resistance to pollination networks, acting as a buffer against species loss. Beside the increasing trend of plant species selectiveness from the seashore inland, we found a significant negative correlation between values of weighted closeness centrality (wCC), as well as of nestedness contribution (n i ), and the mean log-transformed distance from the sea. Animal-pollinated species growing in the pioneer communities of drift line and shifting dunes were more central in the topography of pollination networks (i.e., keystone species 47 ) than those occurring more inland (i.e., peripheral species 47 ). According to Martin González et al. 15 the removal of keystone plant species would have the strongest effects on a pollination network and might even lead to its collapse by triggering the local secondary extinction of specialized pollinator species.
Coastal dune pollination networks were non-random during the overall flowering season. Among network parameters, a non-random nested organization of pollination interactions is expected to increase system resistance by decreasing the probability of local secondary extinction of specialists, which are considered the most vulnerable network members 59 . Furthermore, a non-random modular organization of pollination interactions has been recognized to slow the spread of disturbances, further increasing the resistance of pollination networks 60 . However, values of network selectiveness, nestedness and modularity revealed that the resistance of coastal dune pollination networks was relatively low when compared with pollination networks assessed in other plant communities 52,61,62 . Such relatively low resistance of coastal dune pollination networks further confirms the crucial role of plant species of drift line and shifting dune communities for the prevention of pollination network collapse.
Observed non-random sea-inland pattern of pollination interactions might be the result of the selective pressure exerted by the environmental disturbance on animal-pollinated species. Indeed, according to the specialization-asymmetry-disturbance hypothesis, both specialization and asymmetry of interactions (i.e., nested pattern of interactions 63 ) would explain species' responses to disturbance 64 . In support to this theory, we showed that plant species selectiveness (d′) increased along the sea-inland gradient, while their nestedness contribution (n i ) decreased. In fact, the environmental disturbance to which drift line and shifting dune communities are subjected may increase the stochasticity of pollination interactions, mostly affecting specialized, highly selective, species 6 . Moreover, we proved that animal-pollinated species growing in proximity of the seashore flowered for a longer period than those occurring in the more inland plant communities. A longer flowering period may, in fact, increase the chances of pollination also under variable environmental conditions and, in turn, assure stability to pollination networks during the overall flowering season 65 .
In conclusion, despite some criticalities, our results provided new insights into pollination interactions along environmental gradients. The decreased completeness of pollination interactions from the seashore inland should be further investigated, to disentangle the effect of undersampling from the spatial pattern of species selectiveness. Furthermore, a trait-based approach, aimed at revealing the sea-inland gradient of floral resources, may improve our understanding of the observed spatial pattern of pollinator species richness.
Nevertheless, our study proved that, when different plant communities are spatially interconnected, the assessment of pollination interactions at the landscape scale results in a better understanding of ecosystem dynamics and resistance to perturbations. In the case of coastal dune systems, our results highlighted the essential function of plant species growing nearby the seashore in conferring resistance to pollination networks. Thus, the key function of plant species of drift line and shifting dune communities is not limited only to the creation and maintenance of dune ridges 19 , but also to the maintenance of the stability of pollination interactions over the flowering season.
In order to increase the effectiveness of conservation practices in coastal dune systems, it will be crucial to consider the structure of pollination interactions at the landscape scale along the sea-inland gradient. Focusing on the key role of plant species of the drift line and shifting dune communities would assure resistance, and thus stability, to coastal dune pollination networks, ultimately increasing the long-term sustainability of conservation actions.

Data Availability
Data analysed during this study are included in the Electronic Supplementary Material 2.