Patterns of variation in plant diversity vary over different spatial levels in seasonal coastal wetlands

To quantify the responses of alpha and beta diversity to multivariate gradients, incorporating variation in environmental and management variability in coastal dune slacks.


| INTRODUC TI ON
Biodiversity is declining at an unprecedented rate (Díaz et al., 2019).
This will have significant global impacts on ecosystem services, which are essential for human health and wellbeing (Dasgupta, 2021). To understand these changes, it requires determining which factors control biodiversity at different spatial, temporal and taxonomic scales (Chase et al., 2019;McGill et al., 2015;Sax & Gaines, 2003).
The processes which control diversity-including extinction, extirpation (regional-level extinction), speciation and immigration-can operate at different scales. At a regional scale, climate, dispersal and biogeographic history can strongly control plant community diversity; at a local scale, species interactions, management and resources can exert stronger control (Harrison & Cornell, 2008;Lortie et al., 2004). These impacts are further complicated by species relatedness (Cavender- Bares et al., 2006). In general, more closely related species respond in more similar ways than less-related species (Burns & Strauss, 2011). Therefore, the responses of biodiversity metrics to environmental gradients at higher taxonomic resolutions (e.g. genus) may have high similarity to those at the species level.
As such, exploring biodiversity at a range of spatial and taxonomic scales can provide key insights for understanding the drivers of diversity. This has important consequences for managing biodiversity loss, as different anthropogenic drivers will impact biodiversity at different scales (Vellend et al., 2017).
The importance of scale is difficult to disentangle because of the practical difficulties in obtaining comparative data at multiple scales (Chase et al., 2018). While there is a large body of literature highlighting the role that spatial scale can play in shaping biodiversity patterns, they tend to focus on a single driver of diversity change.
This provides excellent understanding of how biodiversity is impacted by individual drivers such as along temperature, precipitation or latitudinal gradients (Fitzpatrick et al., 2013;Franklin et al., 2019;Ibanez et al., 2018). However, environmental gradients are multivariate; multiple drivers of biodiversity change operate together. If only individual drivers are explored, interactions between drivers of biodiversity change could be missed. Therefore, understanding of real-world patterns requires consideration of real-world gradients. This is essential if we are to scale up and/or down anthropogenic impacts on biodiversity.
Coastal sand dunes are one of the most at risk habitats in Europe and are important reservoirs of biodiversity (Janssen et al., 2017).
Dune slacks are seasonal coastal wetlands, which in the UK support several species of conservation importance (JNCC, 2021). They are expected to be particularly sensitive to environmental changes due to complex interactions between hydrology, nutrient inputs and plant communities (Curreli et al., 2013;Rhymes et al., 2016). Depth to the water table, seasonal fluctuations and duration of drought or flood can alter species interactions and therefore community assembly Willis et al., 1959). The hydrological regime also influences nutrient levels within the dunes, with high water levels maintaining low nitrogen and phosphorus (Lammerts & Grootjans, 1997). Successional processes are also important in structuring dune slack plant communities. Pioneer species establish on bare soil containing low or no nutrients. Over time, plants adapted to these low nutrients colonize, and a moss layer develops. In mature phases, the soil organic matter layer has developed and nitrogen mineralization increases. There is a decline of typical dune slack plants and increase in taller species as competition shifts from nutrients to light (Berendse et al., 1998;Grootjans et al., 2008;Lammerts et al., 1999). Soil organic matter can alter plant responses, as thicker organic matter layers (indicating older dune slacks) have a higher water holding capacity (Bordoloi et al., 2019). Management, such as grazing, may help to slow down successional processes and therefore alter the plant community composition (Millett & Edmondson, 2013). These complex interactions are disrupted by anthropogenic pressures such as climate change, with species that are already at risk being more likely to decline (Bartholomeus et al., 2011). Nitrogen deposition can also disrupt these interactions, with species richness being lowest in sites receiving high levels of nitrogen deposition (Field et al., 2014).
How each wetland will respond will depend on the geographic location , vertical elevation (Moeslund et al., 2013), hydrological regimes (Silvertown et al., 2015), successional processes Grootjans et al., 1991;Jones et al., 2008) and management . To fully understand these drivers of biodiversity variability requires that we disentangle impacts on both alpha and beta diversity. Drivers of alpha and beta diversity are related because beta diversity results from changes in plant communities in space and time, particularly along environmental gradients. However, they are not necessarily tightly coupled because changes in community composition does not always result in changes in alpha diversity. Therefore, to predict how dune slack plant diversity will respond to the various pressures, understanding is needed of how complex environmental interactions drive biodiversity patterns across different spatial and taxonomic levels. Determining how scale of measurement impacts the interpretation of the drivers of biodiversity is essential for a complete understanding of diversity patterns (Anderson, 2018).
Without such understanding, it is difficult to determine the drivers of diversity, the impacts on biodiversity and potential ways to manage these impacts.
In this study, we investigated the spatial-scale dependence of environmental and management drivers of plant community diversity in coastal dune slacks across the UK, and how the measured response depends on the taxonomic resolution considered.
Specifically, the aim was to quantify the responses of alpha and beta diversity to multivariate gradients, incorporating variation in environmental and management variability in UK dune slacks. To do this, three steps were followed: (1) We identified the key environmental gradients between dune slacks; (2) we explored the multivariate gradients associated with patterns of alpha diversity (diversity in an assemblage) across spatial (coastal sand dune system, dune slack and quadrat) and taxonomic (species and genus) levels; and (3) we explored the multivariate gradients driving beta diversity (the variation in composition of assemblages) across spatial and taxonomic levels.
As beta diversity can reflect turnover (replacement) and nestedness (loss), we explored changes in both turnover and nestedness over spatial and taxonomic levels. We also explored the spatial and multivariate environmental drivers of turnover, and how the rate and intensity of turnover differs along these multivariate environmental gradients.

| Coastal sand dune systems
Coastal sand dune systems were selected to represent geographic, climatic and hydrological variability around the UK (Figure 1a). Longterm hydrological recording was needed to fully understand plant responses to previous hydrological conditions (Silvertown et al., 2015), but only a subset of the sites had a history of hydrological recording (using a dipwell). Therefore, to address this gap, we installed dipwells in four sand dune systems to monitor water levels. Specifically, we installed dipwells in Coll Nature Reserve, Sands of Forvie, Tentsmuir and Winterton dunes between June and August 2017. In total, 41 dipwells (one dipwell per dune slack), within 12 coastal sand dune systems across the UK were surveyed.

| Community composition
Between June and September in 2017 and 2018, 164 quadrats across 41 dune slacks within 12 UK coastal sand dune systems were surveyed in the field (Figure 1a,b). In each dune slack, percentage cover of vascular and non-vascular plant species was recorded in, usually, four quadrats (1 m × 1 m) located 4 m north, east, south and west of a central point (the dipwell). This fixed approach around the dipwell ensured that data collection was consistent across the UK.
The location and elevation of each dipwell (41) and each quadrat (164) were recorded using a differential GPS (DGPS) (R6 Trimble RTK-GPS), with vertical accuracy <10 mm (Trimble, 2013). These data were used to calculate quadrat-level hydrological metrics, described below.
Vascular and non-vascular plants were recorded to species level, nomenclature follows Paton (1999), Smith and Smith (2004) and Stace (2019). Vascular plant height was measured at 10 randomly located points within each quadrat. Plant community survey data were used to calculate three measures of alpha diversity: (1) species richness, (2) Shannon's diversity index and (3) Pielou's evenness index. These measures were calculated using the package "vegan" in R (Oksanen et al., 2015;R Core Team, 2022). Beta diversity was calculated using Bray-Curtis dissimilarity at the levels of sand F I G U R E 1 (a) Geographic locations of the 12 UK coastal sand dunes. The blue grids indicate locations based on 10 km presence of dune slacks in the UK (EU habitats directive annex I habitat-H2190) Data downloaded from Eionet central data repository (2020). Map image source: R package "rworldmap" (South, 2011). (b) Schematic diagram of the different spatial levels analysed (from largest to smallest level): Coastal sand dune system (n = 12), dune slack (n = 41) and quadrats (n = 164) dune system, dune slack and quadrat (Figure 1b), using the package "BiodiversityR" (Kindt, 2018;R Core Team, 2022).

| Atmospheric deposition and climate data
Atmospheric deposition and climate data were collated to determine their potential environmental impacts on species diversity in dune slacks. Modelled atmospheric deposition data were downloaded from the Concentration Based Estimated Deposition (CBED) dataset (Smith et al., 2000). This consisted of 5 km × 5 km gridded averages of oxidized nitrogen (NO 2 /NO 3 ), reduced nitrogen (NH 3 /NH 4 ), nonmarine sulphur (SO 2 /SO 4 ) and base cations (Ca and Mg). These were averaged over a 30-year period (1987-2017) because this has been shown to be the best predictor of plant community response to nitrogen (N) deposition (Payne et al., 2019). Climatic variables were extracted from the HadUK-Grid (Hollis et al., 2020). This consisted of 5 km × 5 km gridded daily averages of mean, maximum and minimum annual temperature, and mean annual precipitation averaged over a 10-year period (2007-2017) (Britton et al., 2017).

| Harmonizing the hydrological data
Water table level for each dune slack was measured from a single dipwell, with depth from the ground surface to the water table recorded using data loggers or a water level sensor. The monitoring of water table level varied for each dune slack, with differences in the frequency of data collection (recorded daily, weekly or monthly) and length of data collection (the shortest monitoring period from 2017 to present and the longest from 1974 to present). To address these differences in recording frequency, fill any data gaps and to allow comparisons of water levels between dune slacks for a common monitoring period, there was a need to harmonize the water table data. To do this, water table levels were modelled using the time series analysis software Menyanthes (Bartholomeus et al., 2008;KWR, 2020).
Menyanthes uses daily climate data calibrated against measured water table depth from each dipwell to derive a daily time series model (von Asmuth et al., 2012). Daily climate data required was mean temperature (°C), precipitation (mm) and potential evapotranspiration (PET) (mm). PET was calculated using the Penman-Monteith formula for a reference vegetation (Allen et al., 1994), function "ET. PenmanMonteith" in the R package "Evapotranspiration" (Guo et al., 2016). Data required for PET included maximum, minimum and average temperature (°C), relative humidity (%) precipitation (mm), daily sunshine hours (h) and wind speed (m s −1 ). Climate data were downloaded from coastal weather stations nearest to each sand dune system (MetOffice, 2022). A daily time series model was created for each dipwell for the period of measured water table depth, and the outputs (calibration statistics) were checked. Based on the daily time series models, the hydrological period could be extended by simulating the water table depths for each dipwell. For instance, hydrological data from 2017 to 2021 could be extended to 2014-2021. Menyanthes provides guidelines on how far beyond the core data period it is appropriate to extend the simulation, depending also on the fit of the model. These simulations were particularly important for sand dune systems where the hydrological recording began in 2017, as the simulated water table depths from 2014 to 2018 were used for statistical analysis (see below). Menyanthes was used to model each of the 41 dune slacks (41 dipwells), which reflects local sand dune system hydrology.

| Hydrological metrics
For each dune slack, mean and standard deviation of water  (Curreli et al., 2013), and wet meadows (Gowing et al., 2002(Gowing et al., , 2005. Mean water table level gives a measure of the relative wetness of each plot, providing an understanding of overall prevailing conditions (Richter et al., 1996). At a given point in time, the water table may be higher or lower than this mean, and the extent of this variability will influence plants through drought or waterlogging. Standard deviation of water table level was calculated to give a measure of this variability, and so the extent to which the plant community experiences more extreme conditions. Dune slacks can vary in size, shape and topography . To account for differences in vertical elevation within the dune slack, quadrat-level hydrological metrics were calculated for all 164 quadrats. These metrics were based on the differences in elevation of each quadrat relative to the dipwell recorded using the DPGS.

| Soil characteristics
A soil core (5 cm diameter, 15 cm in depth) was sampled from the south-west corner of each quadrat. The thickness of the soil organic matter was measured in each sample. This was used as a proxy for the age of the dune slack, with older dune slacks having thicker organic matter layers (Grootjans et al., 2008). Organic and mineral horizons were separated. These samples were dried at 40°C, passed through a 2-mm sieve and weighed. Samples were dried at 105°C for 16 h to remove moisture. Loss on Ignition (LOI) at 550°C for 2 h in a furnace was used to provide an estimate of organic matter, and the difference between LOI at 550 and 950°C (for 4 h) to estimate carbonate content (LOC), respectively. Presented are the results for LOI and LOC for the organic horizon.

| Plant nutrient status (vascular leaf sample analysis)
Leaf samples of common dune slack species were collected. To keep the results as consistent as possible, species from four functional groups (grass, sedge, herb or woody) were sampled. Groups sampled included grasses (Agrostis stolonifera, Holcus lanatus), sedges (Carex arenaria, C. flacca, C. nigra, C. oederi and C. rostrata), herbs (Galium palustre, G. verum and Hydrocotyle vulgaris), and a woody perennial (Salix repens). The species included in these samples varied because the same species were not present in every dune slack, and therefore, a random sample of species present in each group were taken.
Plant samples were taken from the quadrats in each dune slack (total of 164 quadrats) for determination of nitrogen (N), phosphorous (P) and potassium (K) content. These values were used to provide a measure of plant nutrient status at each dune slack. Only leaves with no apparent damage were sampled. Leaves were dried in an oven at 60°C for 3 days and ground to a fine powder prior to analysis by an elemental analyser (Forest Research, 2021).

| Management
Management regimes can vary between coastal sand dune systems and may alter plant community composition (Millett & Edmondson, 2013). Vegetation height can be used as a proxy for the grazing regime in sand dune systems. For instance, shorter vegetation can indicate the presence of higher densities of grazing animals (Pakeman, 2004). To give an indication of grazing pressures for each quadrat, we used the average vegetation height (see earlier section), assuming that high grazing pressure will result in shorter vegetation.
Four quadrats (across multiple sand dune systems) were removed from analysis due to missing soil characteristic data, resulting in 160 quadrats. For vascular plants, a total of 171 species were recorded, but nine species (across multiple sand dune systems) were removed from analysis as they were only identified to genus. As a result, in total 162 species distributed across 106 genera were analysed.

| Environmental gradients
To assess the co-linearity between the 16 environmental variables, we created a correlation matrix using the function "corr" in the package "boot" (Canty & Ripley, 2020). We visualized the data using the function "corrplot" from the package "corrplot" (Wei et al., 2017) (see Appendix S1). To account for these correlated variables and to reduce the dimensionality of the environmental data, we performed a Principal Components Analysis (PCA) of the 16 environmental variables using the function "prcomp" in the package "MASS" (Ripley et al., 2013). To ensure that the environmental variables had equal effect on the ordination, the data was standardized, so it has a mean of 0 and variance if 1 (scale = TRUE) using the function "prcomp" in the package "MASS" (Ripley et al., 2013).
A scree plot was used to determine how many principal components should be used (Holland, 2019(Holland, , 2021. The first four PCA axes were chosen as they explained over 74% of the total variation in the data. There was a substantial drop in variation explained with the fifth axis. Loadings were calculated to understand how the environmental variables contribute to each principal component. A positive loading means the variable is positively correlated with the principal component, and a negative loading means the variable is negatively correlated with the principal component (Holland, 2019). As the principal components analysis was performed in a correlation matrix (Holland, 2019), we set the cut-off for associating variables to axes as the critical value for significance based on a two-tailed test using the Pearson's correlation coefficient (df = 39, p = .05, r = .308). Based on this approach, the cut-off value used was 0.308.

| Alpha diversity
To identify the multivariate environmental gradients (PCA axes) associated with patterns in alpha diversity at different spatial and taxonomic levels, we used linear models and linear mixed effects models. For sand dune systems, a linear model was used to predict alpha diversity as a function of the four PCA axes. A linear model was selected as there was no nesting at sand dune system level.
Separate linear models were run for richness, Shannon diversity and Pielou's evenness for species and genus (six models in total). The full model was run using the function "lm" in the base package "stats" (R Core Team, 2022). To identify the most parsimonious model by AIC, backwards elimination was carried out using the function "step" in the package "lmerTest" (Kuznetsova et al., 2017).
Dune slacks and quadrats were analysed using linear mixed effects models (LMM) as they can partition the total variation into different levels, as quadrats were nested within dune slack, and dune slack within sand dune system. At a (i) dune slack level, sand dune system was included as a random effect and the four PCA axes were fixed effects. At the (ii) quadrat level, dune slack nested within sand dune system were included as random effects and the four PCA axes were fixed effects. The response variable were alpha diversity metrics (richness, Shannon diversity or Pielou's evenness), and therefore, separate models were run for each combination of spatial (dune slack and quadrat) and taxonomic level (species and genus) (12 combinations in total). We used a LMM using the function "lmer" in package "lme4," with a Gaussian distribution and the identity link function for (i) dune slack and (ii) quadrat (Bates et al., 2015). For each combination we compared four models, (a) a single level null model, (b) multilevel null model, (c) random intercept model (intercept of PCA axes to vary randomly with (i) dune slack or (ii) quadrat)), and (d) random slope and intercept model (intercept and slope of PCA axes to vary randomly with (i) dune slack or (ii) quadrat)). Models including a (d) random slope and intercept failed to converge for LMM, so results for only (c) random intercept models are presented. To identify the most parsimonious model, backwards elimination using maximum likelihood (REML = FALSE) was carried out on non-significant fixed effects terms using the function "step" (with reduce. random = FALSE) in the package "lmerTest" (Kuznetsova et al., 2017). The model residuals were checked for overdispersion and heteroscedasticity. The final models were fitted with residual maximum likelihood (REML = TRUE) and presented in

| Beta diversity (turnover and nestedness)
Beta diversity (the variation in composition of assemblages) may reflect turnover (replacement) and nestedness (loss). To identify how much of the plant composition variation is due to turnover or nestedness, and how it differs at spatial and taxonomic levels, we used the package "betapart" (Baselga & Orme, 2012). We partitioned the data into different spatial levels (mean of each sand dune system [n = 12], mean of each dune slack [n = 41] and quadrat level [n = 160]) and into taxonomic level (species [n = 162] and genus [n = 106]). For every combination of spatial and taxonomic level, we used the function "betapart.core" in the package "betapart" to create multiple-site beta diversity measures and pairwise dissimilarity matrices based on the composition data. We then used the "beta.
multi.abund" with the family Bray-Curtis in the package "betapart" to partition this into nested and turnover components (Baselga & Orme, 2012;, which was presented as percentage of turnover. To understand the spatial and environmental drivers of turnover, we used a generalized dissimilarity model (GDM) (Ferrier et al., 2007). GDM is a statistical method that analyses and predicts the spatial patterns of turnover across regions. It models the spatial turnover of composition between all samples as a function of the environmental differences between the locations (Ferrier et al., 2007).
We used geographical distances (Latitude and Longitude) and environmental gradients (PCA axes) as predictors, and Bray-Curtis community distance matrix as response variables (Mokany et al., 2022b;Woolley et al., 2017). To prepare the data for GDM analysis, the "formatsitepair function" in the package "gdm" was used (Ferrier et al., 2007;Mokany et al., 2022a). To calculate the proportion of variance explained and the relative importance of each predictor by climate and geographic distance, two "formatesitepair" tables were created (1)

TA B L E 1 Environmental variables used in the principal components analysis (PCA)
community distance (response), and (2) environmental gradients (predictor) and Bray-Curtis community distance (response). Next, the function "gdm" in the package "gdm" (Ferrier et al., 2007;Mokany et al., 2022a), was used to compute three models: (model 1) environmental and geographic distances as predictors, (model 2) only environmental conditions as predictors and (model 3) only geographic distances as predictors (Legendre, 2008). The proportion of variance explained and the relative importance of each predictor of the three models was then compared following .
This process was repeated for all combinations of spatial and taxonomic levels.
To understand how the rate and intensity of turnover differs along environmental gradients for spatial and taxonomic levels, we extracted the I-spline values from the GDM (GDM outputs from model 1) using the function "isplineExtract" in the "gdm" package (Ferrier et al., 2007;Fitzpatrick et al., 2013;Mokany et al., 2022a).
This process was repeated for all combinations of spatial and taxonomic levels.
At a dune slack level and at quadrat level, species and genus richness, Shannon species and genus diversity, and Pielou's species evenness was lowest in older, wetter dune slacks, with taller vegetation and non-marine sulphur deposition (PC3) (Table 3b,c and Figure 2). Species richness was highest when there was higher rainfall, higher base cation deposition and higher organic matter content (PC2) ( Table 3b,c and Figure 2). However, changes in species richness along PC2 were not significant but were retained in the final model as this produced the best model fit (Table 3b,c).

| Beta diversity (turnover and nestedness)
Essentially, all the variance in the composition dissimilarity matrix for sand dune system, dune slack and quadrat was due to turnover (replacement) rather than nestedness (species or genus loss) ( Table 4).
The proportion of community composition variance explained by turnover was lowest at dune system level, and highest at quadrat level, and was consistently lower for genus than species (Figure 3).
Between sand dune systems the total proportion of variance explaining turnover (environment, geographic distance, and shared environment and geographic distance) was 49% for species and 35% for genus; between dune slacks, the proportion of variance explaining turnover decreased to 27% for species and 20% for genus; between quadrats, the proportion of variance explaining turnover decreased further to 20% for species and 20% for genus ( Figure 3). Additionally, across spatial and taxonomic level, there was a strong spatial structure (proportion of variance explained by shared environment and geographic distance, and geographic distance) (Figure 3).
The relative importance of the impact of each environmental gradient on turnover differed between spatial and taxonomic level ( Figure 3). At the sand dune system level, PC4 explained the highest proportion of species and genus turnover. At a dune slack level and at quadrat level, PC3 explained the highest proportion of species and genus turnover ( Figure 3).
The rate and intensity of turnover differed with spatial and taxonomic level, being strongest at quadrat and species level (Figure 4).
Across all levels, the impact of distance between quadrat/dune slack/sand dune system on turnover plateaus is around 400 km. This indicates that with increasing distances up to 400 km the identity of species/genus in a dune slack changes, but over 400 km there is no further change. Across all spatial levels, the amount of composition turnover and rate of composition turnover was low for PC1 and PC2 ( Figure 4). The most important variables in driving plant species and genus dissimilarities at quadrat and dune slack were PC3 (Figure 4).
The highest variable in driving plant species and genus dissimilarities at sand dune system were PC4 (Figure 4).

| DISCUSS ION
Patterns of plant biodiversity in dune slack communities varied depending on the spatial level at which data are collected and analysed, but not the taxonomic level. At smaller spatial levels (dune slack and quadrat), alpha and beta diversity varied along gradients driven by soil characteristics, water table depth and atmospheric deposition. At larger spatial levels (sand dune system), however, patterns of beta diversity were governed by plant nutrient status. Patterns of change in alpha diversity were only detected at some spatial levels, but not others. Indices of diversity (richness and Shannon's) varied at smaller spatial levels (slack and quadrat) while evenness varied mainly at the sand dune system level. The taxonomic level considered made little difference in understanding these patterns.
These results have important consequences for understanding and managing biodiversity because patterns and drivers of biodiversity change can be missed by collecting and analysing data at only one spatial level.
The problem of pattern and scale in biodiversity measurement is well established (Levin, 1992). The environmental drivers of changes in diversity can operate at different scales (Sax & Gaines, 2003) because impacts on the key processes of colonization and extinction are expected to result in complex scaledependant responses (Chase et al., 2019). Thus, the complex scale dependence of biodiversity patterns we have found is not surprising. The consequences of this problem, however, are rarely quantified robustly using real biodiversity data. We found that if surveys were conducted at a sand dune system level (i.e. using a mean community composition per sand dune system) changes in evenness but not richness were detected, whereas if surveys were conducted within a sand dune system (i.e. dune slacks or quadrats) changes in richness but not evenness were detected. Similar differences were measured for beta diversity. At sand dune system level, turnover was driven along a plant nutrient status gradient, whereas within a sand dune system (i.e. dune slacks or quadrats), turnover was driven along an atmospheric deposition, soil characteristics and water table depth gradient. The consequence of this is that management or policy action taken based on data from one spatial level would not necessarily have resulted in optimal outcomes for the biodiversity of the system.
Diversity patterns across two taxonomic levels were explored to understand if species and genus have similar ecological responses.
Phylogenetic correlations between species mean that genus level analyses might provide a fairer test for understanding diversity patterns along environmental gradients (Bonte et al., 2012). Our results show that simplifying the taxonomic resolution resulted in little information loss; the proportion of variation explained by species and genus was similar, largely because most genera were represented by a single species. These findings have an important applied management aspect to dune slacks and demonstrate that identification to a genus level still picks up similar patterns of community turnover as with species. This may save time and money, as in-depth surveys may not be required. However, this generalization may not apply to all habitats as some species-rich habitats contain genera with multiple species within which environmental niches may vary between species in the same genus.
This is the first time that the impact of multivariate environmental and management conditions on patterns of diversity in dune slacks have been explored at different spatial and taxonomic levels. The patterns of variation of plant community composition at smaller spatial levels (dune slack and quadrat) are consistent with some general findings in dune slacks-that diversity is lowest in dune slacks that are older (thicker soil organic matter) and have taller vegetation (lower grazing pressure) (Millett & Edmondson, 2013). Successional processes affect plant community structure. In general, pioneer species establish on bare sand, resulting in initial increases in diversity. Over time soil organic matter and nutrients accumulate, causing a shift towards shrubs and grasses and more intense competition, which results in the extinction of less-competitive species (Grootjans et al., 2008). In older soils with progressive leaching and decalcification, the soils also become more acidic (Berendse et al., 1998;Grootjans et al., 2008), resulting in conditions which support relatively low numbers of species. The grazing regime can also explain these patterns of variation (Millett & Edmondson, 2013 appeared more complex as the hydrological regime influences plant community structure directly though disturbances such as periods of prolonged flooding, and indirectly through mechanisms such as nutrient levels (Grootjans et al., 2008;Jones et al., 2006;Lammerts & Grootjans, 1997). In wetter dune slacks, soil organic matter accumulation is faster (Jones et al., 2008). Along this same environmental gradient of wetter and older dune slack with taller plants, diversity measures were lowest when non-marine sulphur deposition was high. Sulphur emissions have decreased since the 1970s, and may no longer pose a long-term threat to diversity (Smith et al., 2001). Therefore, these results may be in part related to the old legacy gradient of a much higher deposition (Guerrieri et al., 2011;Pakeman et al., 2016).
At higher spatial levels, patterns of beta diversity (turnover) were a consequence of plant nutrient status. Nutrients can alter plant diversity. For instance, in a long-term experiment Storkey et al. (2015) demonstrated that plant diversity increased with decreasing N fertilizer inputs and atmospheric deposition; Ceulemans et al. (2014) demonstrated across European grasslands that plant species richness was negatively related to soil phosphorous. Therefore, it was unsurprising that we observed turnover in dune slack plant communities along a gradient of high to low plant nutrient status. It is also TA B L E 3 (A) Linear model of mean Pielou's evenness at sand dune system level against mean principal components analysis (PCA) axis for sand dune systems. Presented is the final model, after model reduction. important to note that we sampled vascular plant leaves to indicate plant nutrient status within a dune slack, but the species sampled varied between dune slack. This gradient could, therefore, be an artefact of the data collection, reflecting differences in leaf nutrient content between species rather than between dune slacks. However, for each dune slack we sampled from a range of plant species and F I G U R E 2 The relationship between alpha diversity metrics (richness, Shannon diversity and Pielou's evenness) and environmental gradients. Presented are changes in diversity indices along PCA axes 2 (precipitation, atmospheric deposition and soil characteristics), 3 (atmospheric deposition, soil characteristics and water table level) and 4 (plant nutrient status), for different spatial and taxonomic levels.
The spatial levels are sand dune system, dune slack and quadrat, and the taxonomic levels are species and genus. We have presented the separate linear models and linear mixed effects models in one plot for visualization. The points (grey •) represent the 160 quadrats. For sand dune system, the fitted lines (dark blue -----) show the predictions from the linear models. For dune slack (grey _____ ) and quadrat (light blue ----), we have presented the partial plots of the random intercept linear mixed effects models. Changes in species richness along PC2 were not significant but were retained in the final model as this produced the best model fit ensured that this sample covered the same plant functional groups.
We are, therefore, confident that this represents a gradient of plant nutrient status.
Biogeography and the environment can play a key role in species diversity and community structure (Sax & Gaines, 2003). We found that there was a strong spatial signal at sand dune system, dune slack and quadrat level (proportion of variance explained by shared environment and geographic distance, and geographic distance). This suggests that spatially structured environmental factors within a sand dune system and within a dune slack have an important impact on determining community turnover. For instance, within a sand dune system, successional processes, hydrological variability and nutrients can vary considerably (Jones et al., 2006). Within a dune slack, spatial variation in topography means that the water table depth from ground surface varies as vertical elevation varies. Several authors have already noted that fine-scale change in hydrology and nutrient levels structure plant communities (Curreli et al., 2013;Rhymes et al., 2014Rhymes et al., , 2016Silvertown et al., 2015). Community turnover may be also affected by the dispersal ability of different plant species (Bonte et al., 2012), and therefore, the priority effects (the order or the timing of species arrival may affect the growth, establishment or reproduction of a later arriving species) may also impact community structure (Weidlich et al., 2021

F I G U R E 3
The taxonomic turnover along PCA environmental gradients using a generalized dissimilarity model (GDM) of Bray-Curtis distance matrix of abundance data for different taxonomic levels (species and genus) and for different spatial levels (sand dune system, dune slack and quadrat). Sand dune system level beta diversity is calculated from the mean composition of quadrats within the same sand dune system (n = 66), dune slack level beta diversity from the mean of quadrats located within the same dune slack (n = 820) and quadrat level beta diversity from all pairs of quadrats (n = 12,720). For taxonomic and spatial levels, the left-hand bar of each pair represents the proportion of variance explained by only environment, by only geographic distance, their shared variance and unexplained variance. Spatial structure is the combination of geographic distance (black), and shared environment and geographic distance (dark grey). Environmental variation is the combination of shared environment and geographic distance (dark grey), and the environment (light grey). The righthand bar represents the relative importance of each of the PCA axes of the environmental variables. PC1 reflects the differences along a temperature, atmospheric deposition, and water table gradient. PC2 reflects the differences along a precipitation, atmospheric deposition, and soil characteristics gradient. PC3 is gradient of atmospheric deposition, soil characteristics and water table depth. PC4 reflects the differences along a plant nutrient status gradient level (Bossuyt et al., 2003). Future studies should identify the extent to which this missed variation in the environmental conditions experienced by plants is important for understanding patterns of biodiversity.

| CON CLUS ION
Our research demonstrates that diversity patterns in dune slack plant communities vary depending on the spatial level considered, but that the taxonomic resolution considered made little difference in interpretation. Therefore, monitoring of biodiversity change and its drivers must consider variation due to the spatial level measured. If we do not consider patterns across different spatial levels, key environmental and management drivers could be missed. This could have important consequences for the conservation of these European protected habitats which are in decline. In addition, such scale dependence needs to be explored for other measures of biodiversity, such as genetic diversity. Overall, in the effort to understand global biodiversity loss, the importance of scale should take a higher priority than it currently does.

F I G U R E 4
The taxonomic turnover (of species and genus) along PCA environmental gradients using generalized dissimilarity model (GDM) for spatial level (sand dune system, dune slack and quadrat). GDM models-fitted I-splines of Bray-Curtis distance matrix of abundance data for variables that are significantly associated with sand dune system, dune slack and quadrat. The height of the curve indicates the amount of composition turnover and the importance of the axis in explaining beta diversity. The shape of the curve indicates the rate of composition turnover along the PCA axis gradient. The length of the gradient depends on the data; that is, sand dune system is a mean of the 12 sand dunes. Geographic distance reflects the impact of distance. PC1 reflects the differences along a temperature, atmospheric deposition, and water table gradient. PC2 reflects the differences along a precipitation, atmospheric deposition, and soil characteristics gradient. PC3 is gradient of atmospheric deposition, soil characteristics and water table depth. PC4 reflects the differences along a plant nutrient status gradient