Can Vegetation Removal Successfully Restore Coastal Dune Biodiversity?

: Coastal dune habitats have been declining globally over the last several decades due to rapid urbanization. Within remaining dune systems, dune ﬁxation has resulted in further losses of mobile dunes with negative impacts on their associated species. Some studies suggest vegetation removal can initially promote habitat heterogeneity, and increase availability of suitable habitats for psammophile, xeric and endemic mobile dune species, but longer-term responses are generally unknown. We investigated the temporal trends of four taxonomic groups to determine the e ﬀ ect of vegetation removal on dune assemblages over a 12-year period at an LTER site. Three di ﬀ erent forms of removal are investigated here—removal in a grid form on ﬁxed dunes, removal of the wind-facing slope vegetation on semi-ﬁxed dunes and opportunistic o ﬀ -road driving on disturbed dunes. Results were varied across taxa, highlighting the need for multi-taxa monitoring in conservation and restoration management. Overall, ﬁxed dune treatment had very little e ﬀ ect, while a stronger response was found in semi-ﬁxed treatments in particular for mobile dune indicator species, which showed evidence of recolonization within a few years following treatment. Disturbed dunes were most similar to mobile dunes for animal taxa indicating that pulse removal may not be as e ﬀ ective as continuous press disturbance. Nevertheless, a less destructive form of disturbance such as re-introduction of grazing might be preferable and requires further investigation.


Introduction
In the last 70 years there have been extensive losses of coastal habitats and naturally functioning coastal dune systems globally [1][2][3][4][5][6][7]. This is due to dramatic population expansion and urban development along coastal plains, coupled with changes in land-use practices and climate-related changes in aeolian processes. [1,[8][9][10]. Globally this has led to substantial increases in dune fixation (stabilization) and a significant loss of naturally mobile (shifting) dunes [11][12][13][14][15][16]. Restoration strategies for coastal dune systems are complex because sandy dunes are exposed to multiple dynamic processes that operate on different temporal and spatial scales [17,18]. Moreover, dune restoration methods are varied and sometimes conflicting. These range from re-establishment of vegetation, to removal of vegetation and eliminating exotics; the strategy usually depends on the management objectives (reviewed by [2,19]). Many dune management regimes around the world have focused on promoting dune fixation through the introduction of (sometimes non-native) grasses and nourishment in order biodiversity and improve ecosystem functioning, the recovery of biotic components of ecosystem processes across multiple trophic levels must be understood [63][64][65][66]. Restoration of coastal dune systems should therefore consider the response of both faunal and floral species.
Biotic responses are often not reported in landscape restoration, and rarely for animal taxa [59,67]. Most coastal restoration projects only consider geomorphological and vegetation components, and usually are only monitored over a short time frame. Only a handful of coastal restoration studies evaluate more than one taxon concurrently, and a recent review (Lithgow et al. 2013) found that 88% of coastal restoration studies only considered plant responses, usually only in terms of vegetation cover or species richness, not composition. There is a clear need for a coastal restoration studies involving vegetation removal to consider responses of species composition over a long timeframe and across multiple taxa (particularly faunal taxa). We found only one study that investigated composition responses of multiple taxa (arthropods and plants) to disturbance in coastal dunes [49], while Kutiel et al. [68] are the only authors to consider animal and plant responses to vegetation removal concurrently. The Nizzanim Dune Nature Reserve (NDNR), on the southern Mediterranean coast of Israel, provides a rare opportunity to conduct a restoration study considering all these components.
Livestock grazing and firewood collection by nomadic tribes was prevalent for hundreds of years in what is now NDNR [69][70][71]. Exclusion of this anthropogenic presence since the middle of the last century has led to substantial shrub encroachment and dune fixation across the reserve [69,72]. It is expected that in the absence of disturbance, the dunes in NDNR will eventually homogenize and become fixed, which will result in the loss of mobile dunes and their associated species, and an overall loss of βand γ-diversity across the reserve [26].
To mitigate the threat of dune fixation in NDNR, the Nizzanim Long Term Ecological Research (LTER) project implemented two restoration experiments involving perennial vegetation removal as a form of pulse disturbance. The goal of these trials was to restore the heterogeneity of dune states and increase γ-diversity by providing a suitable habitat for mobile dune species [69]. This is the largest landscape scale restoration program undertaken in Israel, and the long-term responses have been monitored as part of the LTER program since 2005.
Parallel to the trials, off-road vehicles (ORV) have been utilizing several dunes in the reserve providing an opportunity to examine the effect of a continuous press disturbance on species composition. In this study, we investigate the effect the two controlled pulse removals and the un-controlled ORV press disturbance, in order to determine whether these different types of vegetation removal recreate suitable habitat for mobile dune species. Restoration can take several years, even decades, and can rarely precisely replicate an original status [2,73]. Therefore, it was expected that by removing vegetation and creating disturbance, the dunes would become more suitable for mobile-dune obligate species, and less suitable for fixed-dune species.
Employing the terminology of Lake [73], we used a M-RCI experimental design; comparing multiple repeated replicates in Reference (mobile), Control (un-treated) and Intervention (treated) dunes.
We addressed the following questions: i. How do species compositions of different taxa respond to different removal treatments? ii. Can indicator species represent assemblage level responses? iii. Can the different forms of removal contribute to successful coastal dune management?

Hypotheses Tested
First, we tested the hypothesis that in control and reference dunes, fixed dunes would support a greater abundance (H 1 ) and species diversity (H 2 ) under productivity theory [74].
Next, we tested the responses to treatment predicting that: H 3 : Reduction of primary productivity through mechanical removal will shift a) abundance, and b) diversity towards reference dune levels Appl. Sci. 2020, 10, 2310 4 of 33 H 4 : Composition of treated dunes will shift away from the control dune composition towards reference (mobile) dune composition for each taxon Over time, responses in composition could therefore be expected to either: H 4 a: Gradually become more similar to the reference dune state (turnover) H 4 b: Present an initial shift and then gradually return to control state (resilience) or H 4 c: Remain the same and not shift away from the control state (resistance) Natural coastal dunes are dynamic habitats exposed to a variety of press (wind erosion/sand burial, grazing) disturbances, which creates assemblages of highly adapted mobile-dune species. Therefore, we predicted that: H 5 : Mechanical press disturbance (ORV disturbance) will be more effective than a single mechanical pulse disturbance (removal by bulldozer), at creating community composition similar to the reference dune.
Finally, we considered the impacts of treatment on individual Indicator Species (IS). Since IS will tend to be highly affiliated to their specific habitat type, species abundance should increase with habitat availability.
H 6 : Treatment of semi-fixed and fixed dunes will result in an increase in abundance of species that are affiliated to mobile dune habitat, while species that prefer semi-fixed and fixed habitats will decrease.

Study Site
The Nizzanim LTER is an on-going collaborative project established in 2004 to monitor plants, arthropods, small mammals and reptiles [26]. The study was conducted at the Nizzanim LTER site in Nizzanim Dunes Nature Reserve (NDNR), Israel (31 • 42'-31 • 44'N, 34 • 35'-34 • 36'E) covering an area of 20 km 2 along the Mediterranean coastline. Annual average temperature is 20 • C, and rainfall is 400-500 mm per annum, falling mainly during winter (November-April). The common wind direction is south-west with a very low Drift Potential Index [70].
Sand dunes in NDNR are a globally distinct ecosystem situated at the intersection between desert and Mediterranean coastal dune systems [69]. The edaphic conditions and the physiognomy of the sand dunes enable numerous Saharan taxa to range northward in spite of the temperate climate that favors Mediterranean species [75], creating an unusual and highly diverse assemblage compared to other Mediterranean coastal systems [76]. Three distinct fixation states of dunes are recognized in NDNR, each state with its own characteristic species diversity and composition [26,77]. Following the declaration of the reserve, livestock grazing and wood harvesting by local communities (predominantly Bedouin) were prohibited and rapid dune fixation occurred, increasing the shrub cover from 4% to 20% of the landscape over a period of 60 years [27,72,78].
The Nizzanim LTER site consists of mobile, semi-fixed, and fixed dunes separated by inter-dune depressions [79]. The three fixation states can be categorized based on the gradient of perennial percentage vegetation cover (PPC), aeolian sand movement, and visual indicators such as dune geomorphic structure, plant distribution, dominant perennial species and soil color [27,56,[79][80][81].
Mobile dunes have 5-15% PPC, mainly Ammophila arenaria, and make up approximately 20% of the reserve [27]. Wind-driven shifting sands result in frequent burial and exposure of perennial plants on mobile dunes. Semi-fixed dunes have 16-30% PPC, are dominated by wormwood (Artemesia monosperma) and desert broom (Retama raetam), and cover approximately 70% of the area. Lastly, fixed dunes vary between 30-50% PPC, and compose approximately 10% of the area of the reserve.

Removal Experiments
Three separate treatments were tested ( Figure 1). For each set of treated dunes, there were paired untreated control dunes and mobile dunes were considered as the restoration reference for all treatments (Figure 1e) (other dunes were monitored under the LTER but were not included in this study). A summary table of the overall experimental design has been included in the supplementary information (Table S1). design was chosen in order to leave some vegetation needed for animal survival, while more closely emulating the distribution of natural vegetation found on mobile dunes [see 77]. Four untreated semi-fixed dunes were selected as controls (Figure 1ii). iii.
Disturbed Dunes Treatment: Off-road vehicle (ORV) disturbance has been occurring illegally in the reserve for many years (Figure 1vi), which began prior to the initiation of the experimental plots. We began monitoring these dunes in 2012, in order to study the impact of ORVs on dune biodiversity. We consider these dunes as 'treated' with an unquantified press disturbance. The exact natural state of these dunes prior to the disturbance is unknown; on the one hand, one might assume they were similar to mobile dunes based on historical aerial photos from 1965 and 1999 [69]. On the other hand, most dunes in the immediate vicinity are semi-fixed. Disturbed dunes were therefore compared with both mobile and un-treated semi-fixed dunes as the reference dunes. Annual plants are almost entirely absent from disturbed dunes and as such this taxon was not monitored on these dunes. Mobile dune (reference) dune; (f) Off-road vehicle disturbance. Semi-fixed and Fixed Treated dunes were paired with un-treated Semi-fixed and Fixed untreated control dunes respectively. All treated dunes were compared to their un-treated control pair and Mobile dunes as the restoration reference. Disturbed dunes were exposed to illegal, un-quantified off-road vehicle disturbance. They were compared to both Mobile and Semi-fixed untreated dunes to examine their relative composition. Frames of the different dune types shown are colour consistent throughout all figures.
i. Fixed Dune Treatment: In October 2005, perennial vegetation was removed from the entire fixed dune in a grid formation using a bulldozer, to accurately reach the desired amount of remaining vegetation (Figure 1b). The removal reduced the vegetation from 31-50% to approximately 15-20% for each dune [69]. Three dunes were treated and three control dunes were monitored in parallel. ii. Semi-Fixed Treatment: In November 2012, an experiment was conducted on four semi-fixed dunes, which had a range of 16-30% average PPC prior to treatment. This time, the experiment involved removal of all vegetation from the wind-facing slope of four semi-fixed treated dunes. The vegetation remaining on the crest and slip-face was left intact (Figure 1d). This experimental design was chosen in order to leave some vegetation needed for animal survival, while more closely emulating the distribution of natural vegetation found on mobile dunes [see 77]. Four untreated semi-fixed dunes were selected as controls (Figure 1c). iii. Disturbed Dunes Treatment: Off-road vehicle (ORV) disturbance has been occurring illegally in the reserve for many years (Figure 1f), which began prior to the initiation of the experimental plots. We began monitoring these dunes in 2012, in order to study the impact of ORVs on dune biodiversity. We consider these dunes as 'treated' with an unquantified press disturbance. The exact natural state of these dunes prior to the disturbance is unknown; on the one hand, one might assume they were similar to mobile dunes based on historical aerial photos from 1965 and 1999 [69]. On the other hand, most dunes in the immediate vicinity are semi-fixed. Disturbed dunes were therefore compared with both mobile and un-treated semi-fixed dunes as the reference dunes. Annual plants are almost entirely absent from disturbed dunes and as such this taxon was not monitored on these dunes.

Rodents
Rodents were monitored in autumn (September-October) from 2005-2016, using 27 Sherman traps per dune in three rows of nine traps each, placed alternately in open and bush patches spaced approximately 10 m apart. Trapping took place over two consecutive nights, using mark recapture to calculate the estimated abundance for each species. We used the Chapman estimator, which is an adjusted form of the Lincoln-Petersen Index [82].

Reptiles
Reptiles were sampled in autumn 2005-2016 using a combination of several methods, including ten pitfall traps (10 L buckets), two track transects (100 m long), four activity transects (100 m long) and opportunistic sightings. Ranked abundance for each species was then calculated using a max-log algorithm across all sampling methods with a range of 0 (absent) to 5 (highly abundant). For full calculation methods see Shacham and Bouskila's work [28].

Beetles
Ground-dwelling beetles were monitored in spring (March-May), between the years 2006-2016 (not all years were sampled). Beetles were collected using dry pitfall traps of 12 cm depth and 10 cm diameter, in a regularized pattern of 10 traps on each dune section (windward slope, crest and slip-face), alternating between open and shrub patches. Pitfall traps remained open for 36 hours. The 30 traps were then pooled to give one dune sample.
Beetles were trapped live, identified to morphospecies (sensu [83]), and released where possible or collected and later identified by experts. Some species remained as morphospecies or were identified to genus or family, i.e., to the best Recognizable Taxonomic Unit (RTU) possible. Since data was collected using only one method, we were able to use a direct measure of abundance for beetles (compared to the rank abundance used with the mixed methods for reptiles, or estimated abundance for mark-recaptured rodents).

Annual Plants
Annuals were monitored in spring (March-May) 2006-2016, using 40 × 40cm quadrats. As for beetles, all three sections of the dune were sampled, resulting in a total of 72 quadrats per dune. In each slope a 100 m 2 plot was marked out, and within each plot 12 open and 12 bush quadrats were randomly placed. Total dune cover was then standardized based on the area contribution of each slope to the whole dune.

Data Analysis
For all analyses we tested each taxon separately. First, species with a total abundance of <5 individuals (or 5% cover for annuals) across all years and all dunes were removed. We also removed all species that had been found in less than three years, irrespective of their total abundance. Each experimental treatment type was then examined separately; for each we compared the treated dunes to the respective control dunes (e.g., treated and un-treated fixed dunes) and the reference dunes (mobile). In a similar manner, for the disturbed dunes we consider un-treated semi-fixed dunes as the control pair, and mobile dunes as the reference. All analyses were performed in R version 3.4.3 [84].

Assemblage Abundance and Richness
To understand the effects of treatment on assemblages, we examined the response of assemblage abundance, species richness and composition. First, we compared control dunes (untreated semi-fixed and fixed dunes) with the restoration reference (mobile dunes). Then we used these differences to set the expectations for restoration effects. For example, if there were no differences between control and reference states, we expected no difference between treatment and control, or treatment and reference states (effectively H 1 would be the same as H o ). Conversely, if there were differences between control and reference, then treated dunes were expected to become significantly different from the control, and would either become more similar to the reference, or different from both control and reference.
For each treatment we tested the treated dunes against their untreated control and against the mobile dune reference. Estimated species richness, abundance of rodents, and percentage cover of annuals are all continuous, numerical variables. Therefore, to test differences in these parameters we used linear mixed-effect models as described in Laird and Ware [85]. Since dune sampling was not balanced across all years, we set year as a random intercept effect without interactions for any models, using the nlme package [86].
Generalized Linear Mixed Models (GLMM) allow the introduction of random effects in hierarchical, non-linear, count data [87]. Species richness of all taxa, total abundance of beetles, and rank abundance of reptiles are count data. Therefore, we tested these parameters for differences among dune states using GLMM with Penalized Quasi-Likelihood (GLMM-PQL) with the MASS R package [88], with dune state specified as a fixed effect and year as a random intercept effect, under a Poisson distribution.

Composition
To test the effect of treatments on dune composition we used Redundancy Analysis (RDA) in the R package vegan [89], in a constrained model with year specified as a conditional factor: rda(species matrix)~dune.treat * year + Condition(year) (1) This allowed us to examine the effect of dune type (i.e., treated, reference and control dunes) after the effect of year had been removed as a covariable. We then used the scores along the first RDA axis (RDA1) to calculate the treatment effect as a proportion of total distance between control and reference dunes, shown by directional arrows. ). Finally, we used pairwise Permanova to test the effects of treatment on the centroids of treated dunes compared to their control pairs only (i.e., without mobile dunes), with years as replicates. Multivariate methods were redundant for examining rodents (since there are effectively only two species), but we included ordination plots for consistency in comparison with other taxa.
To visualize the multivariate changes over time, we plotted the RDA ordination as a function of year and used the mean scores for each dune type on the first RDA axis (RDA1). This is conceptually related to the Principal Response Curves approach (PRC) [90], but here our analysis is conducted without standardizing to a baseline. We did this because we wanted to be able to see the fluctuation in all dune types in response to environmental stochasticity.
Since the data were unbalanced, we could not run a standard PRC permutation test, so to test for temporal effects on treated dunes, we used Permutational Multivariate Analysis of Variance (Permanova). These tests were conducted for each year separately, with the treated dune type set as the reference dunes with which to compare both the control pair and the reference (mobile) dunes simultaneously.

Indicator Species
We identified Indicator Species (IS) for each taxon in two stages. First, we identified species that are highly affiliated to a single dune state (mobile, semi-fixed or fixed dunes) using the IndVal method. This method was proposed by Dufrêne and Legendre [91] to identify the ecological preference of a given species among a set of alternative site groups. An Indicator Value (IV) is a measure of association ranging from 0% (individuals are spread equally between the dune types) to 100% (all individuals of a species are observed in all dunes of only one dune state i.e., the indicated dune state). IVs were calculated for each species with the indicspecies package in R [92]. To avoid rare species that were unlikely to be captured or useful for monitoring, we only included common species with total abundance of 20% or more of the average species abundance within each taxon. Note this was in addition to the species that were already removed from the whole study as described earlier in Section 2.4. IVs were tested against a null distribution using a permutation test [93]. IVs that were significantly different from the permuted random distribution (p < 0.01) were considered as strong candidates to be indicator species.
The value of IV is the degree of affinity of each species to a given dune state, but it does not indicate the vegetation cover that the species is associated with. To understand this association, we calculated a second measure which reflects the percentage of vegetation in which a species is most commonly active; this is effectively PPC weighted by abundance (wPPC) such that: where n is the total abundance over time of species i in dune j, and PPC j is the average perennial plant cover in dune j. Thus, a highly obligate mobile dune species would have a high IV and a low wPPC, while a highly obligate fixed dune species would have a high IV and a high wPPC.
Theoretically, wPPC can range from 0% to 100%, but the dunes in the Nizzanim TER varied from 10% to 36% cover, so these are the minimum and maximum possible values in this study, respectively. Equation (2) forms the preliminary steps for calculating the Species Sandiness Index (SSI) in Rubinstein et al. [77], but the range of SSI is dependent on the range of plant cover included in the study and is designed to provide a scale where sandiness in mobile dunes has higher values than in dunes with high vegetation cover. The value of wPPC for each species is comparable across other habitats, and the positive direction reflecting increasing PPC is more informative in the context of this study.
For each taxon, we objectively selected the two IS with the lowest and highest wPPC, to reflect the species most affiliated to mobile and fixed dune states respectively. This stage was done using data from control dunes only. These IS are expected to be the most specialized, adapted or obligate species for the respective dune states. That said, on the fixed dunes they could also be generalist species that are at the edge of their natural (e.g., inland) range and unable to inhabit the more dynamic dune states.
We then tested whether effect of the treatments on the eight selected IS was according to the expectations, i.e., away from the control dunes towards the reference dunes. According to the H 1 hypothesis, we should expect that mobile IS would increase in abundance in the treated dunes, while fixed IS would decrease.
The abundance of each IS was tested using GLMM-PQL (reptiles and beetles) and LME (rodents and annuals) models similar to those used for assemblage abundances [26]. We assessed the temporal trends of each IS using a permuted Empirical Cumulative Distribution Function (ecdf ) in the R stats package [84]. This simple permutation approach generates a p-value for the null hypothesis of zero difference in abundance of a species between two habitats. It allows us to consider unbalanced data with zero-inflated distributions that are not testable under ANOVA. We conducted each test on pairwise treatment and control, or treatment and reference dunes for each year separately.

Abundance & Richness
We tested assemblage-level abundance and species richness to determine whether these parameters were affected by vegetation removal. We considered any result with p < 0.05 as significant Full results and p-values are given in Supplementary Tables S2 and S3. We plotted the results as boxplots where the upper and lower edges of the box correspond to the 25% and 75% quartiles (Figures 2 and 3). The upper whiskers extend from the edge of the box to the highest value, excluding outliers. Data beyond 1.5 of the 25-75% inter-quartile range are outliers and plotted as points. Note that for many of the rodent samples, either one or two species was found. Therefore, for some dune types, the boxplot appears as a single line with no variation. On its own, boxplots for the rodent data would normally not be the most appropriate way to present these results. However, we felt it pertinent to maintain the same graphics across taxa to make them more easily comparable.
Appl. Sci. 2020, 10, x 10 of 34 between treated and untreated fixed dunes, both remaining significantly different from mobile dunes (Figure 3i,l) and contradicting the expectations for a decrease (rodents) and increase (annuals) in response to treatment. The apparent reduction in reptile richness in treated dunes became similar to (i.e., not significantly different from) to mobile dunes (Figure 3j), but the decline was not large enough to become significantly different from the untreated control dunes (p = 0.10) so the expected result was only half met. Beetles showed a small increase in treated fixed dunes towards mobile dunes compared to untreated fixed dunes, and thus beetle species richness did not differ significantly from the reference or the control ( Figure 3k). Lastly, dunes disturbed by ORV had significantly higher rodent species richness compared to untreated semi-fixed dunes, similar to mobile dunes ( Figure 3m). Reptiles were found to have significantly lower species richness on disturbed dunes compared to semi-fixed control dunes, but were not significantly lower than mobile dunes (Figure 3n). Beetle richness in treated semi-fixed dunes remained similar to richness in both untreated semi-fixed dunes and mobile dunes (Figure 3o). Each treated dune type was tested against its untreated pair (control) and mobile (reference) dunes. Disturbed dunes were not monitored for annuals, and untreated semi-fixed dunes were used as the control. Significant differences were tested with generalized mixed models for beetles and reptiles (Poisson), and linear mixed models for rodents and annuals. p < 0.05 for treated dunes compared to * control and + reference dunes.  d) Each treated dune type was tested against its untreated pair (control) and mobile (reference) dunes. Disturbed dunes were not monitored for annuals, and untreated semi-fixed dunes were used as the control. Significant differences were tested with generalized mixed models for beetles and reptiles (Poisson), and linear mixed models for rodents and annuals. p < 0.05 for treated dunes compared to * control and + reference dunes. Each treated dune type was tested against its untreated pair and reference (mobile) dunes. Disturbed dunes were not monitored for annuals and un-treated semi-fixed dunes were used as the control. Significant differences were tested with generalized mixed models for beetles and reptiles (Poisson), and linear mixed models for rodents and annuals. p < 0.05 for treated dunes compared to * control and + reference dunes, respectively.

Composition
To examine the effect of treatments on composition we used constrained multivariate ordination. Firstly, we conducted a redundancy analysis (RDA) of all dune states together, in order to compare the relative composition of all dunes in the study, and to visualize the relative magnitude of the effect across press and pulse treatment types (Figure 4a-d). We then tested the treatment effect in each dune state separately, in order to meet the assumptions for the statistical tests ( Figures 5-7). For all taxa, composition was a good indicator of dune state as expected under hypothesis H2, as seen by the strong clustering and gradient from fixed control to mobile dunes in all taxa along the first RDA axis (RDA1). Explained variation was much higher in RDA1 than RDA2 for all taxa ( Figure 4).
Rodent assemblages in Nizzanim coastal dunes were heavily dominated (98%) by two gerbil species Gerbillus pyramidum and G. andersoni allenbyi. Therefore, ordination formed an arc effect along the first axis ( Figure 4a) as the ratio between the two species shifted from only G. allenbyi in fixed dunes, to only G. pyramidum in disturbed dunes. The disturbed dune centroid was located beyond (further left) the mobile dune's along RDA1 (yellow arrow, Figure 4a), since both species are present in the latter.
Both semi-fixed and fixed dunes showed minor shifts towards mobile dunes (the expected direction for our restoration goals), as a result of treatment for faunal taxa (blue and green arrows in Figure 4a-c), which supported hypothesis H4. A small and insignificant shift in the opposite direction was seen for annuals (green arrow; Figure 4d). Annuals in both semi-fixed control and treated dunes were similar to mobile dune in RDA1 (Figure 4d). ORV disturbed dune composition was virtually identical to mobile dune composition for beetles (Figure 4c), and was similar for reptiles along RDA1, but was distant along RDA2 (yellow arrow; Figure 4b). Overall, press treatment (ORV) appeared to be more similar to reference dunes, compared to pulse treatments, in accordance with hypothesis H5.   Each treated dune type was tested against its untreated pair and reference (mobile) dunes. Disturbed dunes were not monitored for annuals and un-treated semi-fixed dunes were used as the control. Significant differences were tested with generalized mixed models for beetles and reptiles (Poisson), and linear mixed models for rodents and annuals. p < 0.05 for treated dunes compared to * control and + reference dunes, respectively.
We first tested un-treated fixed control and semi-fixed control dunes against the reference mobile dunes, to set up the expectations for treatment effects. There was no difference between control and reference dunes for any of the faunal taxa (Figure 2a-c). This meant we rejected H 1 that abundance of animals increases with increasing primary productivity. For annual plants, increasing perennial cover led to increased annual cover as expected ( Figure 2d).
These findings in control dunes established the expectations for responses to treatment (H 3 ): For annuals, we expected treatment to result in a decrease in cover, while for animal taxa, no difference was found between control and reference dunes, so no difference would be expected in response to treatment.
For the semi-fixed treatment, there was no significant difference in abundance between treated and untreated semi-fixed dunes or mobile dunes for all three faunal taxa (as expected based on the expectations set up in the previous statement ( Figure 2b,e-g). In the fixed dune treatment, rodents maintained similar (not significantly different) abundance to un-treated control and mobile dunes, as expected (Figure 2i), while reptile abundance was significantly lower in treated fixed dunes than their control pairs ( Figure 2j). Conversely, beetle abundance was significantly higher in treated fixed dunes than both control and reference dunes ( Figure 2k). Thus, for reptiles and beetles, treatment did have an effect, despite there being no expectation for any treatment effect. Disturbed dunes exposed to off-road vehicle disturbance were also expected to be more similar to mobile and less similar to untreated semi-fixed control dunes. There was no effect of vehicular disturbance on rodent abundance ( Figure 2m). However, reptile abundance was significantly lower on these dunes compared to both mobile and semi-fixed control dunes (Figure 2n). In contrast, beetle abundance was significantly higher in disturbed dunes than in mobile and semi-fixed control dunes ( Figure 2o).
Cover of annuals increased with increasing dune fixation so, treatment was expected to reduce % cover of annuals, making the treated dunes more similar to mobile dunes. However, on semi-fixed dunes, annual % cover did not significantly decrease in response to treatment, with % cover remaining similar to semi-fixed controls and significantly higher than in mobile dunes ( Figure 2h). In treated fixed dunes, the % cover was significantly higher than both mobile and fixed control dunes, which meant the direction of response was opposite to expectations ( Figure 2l).
As with abundance, we tested the control dunes for differences in species richness (Figure 3a-d). For rodents and beetles there was a gradient of declining species richness with increasing perennial cover (fixation state) in control dunes (Figure 3a,c), which opposes the productivity-diversity theory (H 2 ). Thus, the expectations for treatment would be an increase in richness following removal, if treatment was successfully replicating mobile dune conditions and therefore biodiversity structure. Note that for beetles, the decline in richness was only significant between mobile and fixed dunes, so the semi-fixed treatment was expected to remain similar to both control and reference. Conversely, reptile and annual plant richness increased with increasing dune fixation (Figure 3b,d), supporting the theory that productivity increases diversity. The expectations for these latter taxa were thus that species richness would decrease in response to treatment.
In semi-fixed dunes, no response of richness was seen in reptiles ( Figure 3f) and annuals ( Figure 3h) which meant an absence of expected effect, while beetle richness remained unchanged in response to treatment as expected (Figure 3g). In fixed dunes, rodent and annual richness showed no difference between treated and untreated fixed dunes, both remaining significantly different from mobile dunes (Figure 3i,l) and contradicting the expectations for a decrease (rodents) and increase (annuals) in response to treatment. The apparent reduction in reptile richness in treated dunes became similar to (i.e., not significantly different from) to mobile dunes (Figure 3j), but the decline was not large enough to become significantly different from the untreated control dunes (p = 0.10) so the expected result was only half met. Beetles showed a small increase in treated fixed dunes towards mobile dunes compared to untreated fixed dunes, and thus beetle species richness did not differ significantly from the reference or the control (Figure 3k).
Lastly, dunes disturbed by ORV had significantly higher rodent species richness compared to untreated semi-fixed dunes, similar to mobile dunes (Figure 3m). Reptiles were found to have significantly lower species richness on disturbed dunes compared to semi-fixed control dunes, but were not significantly lower than mobile dunes (Figure 3n). Beetle richness in treated semi-fixed dunes remained similar to richness in both untreated semi-fixed dunes and mobile dunes (Figure 3o).

Composition
To examine the effect of treatments on composition we used constrained multivariate ordination. Firstly, we conducted a redundancy analysis (RDA) of all dune states together, in order to compare the relative composition of all dunes in the study, and to visualize the relative magnitude of the effect across press and pulse treatment types (Figure 4a-d). We then tested the treatment effect in each dune state separately, in order to meet the assumptions for the statistical tests ( Figures 5-7). For all taxa, composition was a good indicator of dune state as expected under hypothesis H 2 , as seen by the strong clustering and gradient from fixed control to mobile dunes in all taxa along the first RDA axis (RDA1). Explained variation was much higher in RDA1 than RDA2 for all taxa (Figure 4).
Rodent assemblages in Nizzanim coastal dunes were heavily dominated (98%) by two gerbil species Gerbillus pyramidum and G. andersoni allenbyi. Therefore, ordination formed an arc effect along the first axis (Figure 4a) as the ratio between the two species shifted from only G. allenbyi in fixed dunes, to only G. pyramidum in disturbed dunes. The disturbed dune centroid was located beyond (further left) the mobile dune's along RDA1 (yellow arrow, Figure 4a), since both species are present in the latter.            Figure 5. Redundancy analysis for (a) rodents, (b) reptiles, (c) beetles and (d) annual plants in treated semi-fixed dunes (pale blue triangle) compared to un-treated semi-fixed (control) dunes (dark blue square) and Mobile dunes (orange square) which are the treatment reference. Distance bars along the 1st RDA axis show the treatment effect (blue arrow) as a percentage of total distance between control and reference (black dashed). Significant differences were tested using pairwise Permanova. * and + represent p < 0.05 for treated dunes compared to control and reference dunes, respectively. and reference (black dashed). Significant differences were tested using pairwise Permanova. * and + represent p < 0.05 for treated dunes compared to control and reference dunes, respectively.

Figure 6. Redundancy analysis for (a) rodents, (b) reptiles, (c) beetles and (d) annual plants in treated
fixed dunes compared to untreated fixed (control) and mobile dunes (reference). Distance bars along the 1st RDA axis show the treatment effect (green arrow) as a percentage of total distance between control and reference (black dashed). Significant differences were tested using pairwise Permanova. * and + represent p < 0.05 for treated dunes compared to control and reference dunes, respectively.   Figure 6. Redundancy analysis for (a) rodents, (b) reptiles, (c) beetles and (d) annual plants in treated fixed dunes compared to untreated fixed (control) and mobile dunes (reference). Distance bars along the 1st RDA axis show the treatment effect (green arrow) as a percentage of total distance between control and reference (black dashed). Significant differences were tested using pairwise Permanova. * and + represent p < 0.05 for treated dunes compared to control and reference dunes, respectively.
Both semi-fixed and fixed dunes showed minor shifts towards mobile dunes (the expected direction for our restoration goals), as a result of treatment for faunal taxa (blue and green arrows in Figure 4a-c), which supported hypothesis H 4 . A small and insignificant shift in the opposite direction was seen for annuals (green arrow; Figure 4d). Annuals in both semi-fixed control and treated dunes were similar to mobile dune in RDA1 (Figure 4d). ORV disturbed dune composition was virtually identical to mobile dune composition for beetles (Figure 4c), and was similar for reptiles along RDA1, but was distant along RDA2 (yellow arrow; Figure 4b). Overall, press treatment (ORV) appeared to be more similar to reference dunes, compared to pulse treatments, in accordance with hypothesis H 5 .
Changes in response to treatment were tested for each dune type separately, comparing each treatment to the reference (mobile) dunes, and their control dunes for semi-fixed ( Figure 5), fixed ( Figure 6) and disturbed dunes (Figure 7). Pairwise Permanova results are given in Table S4. Rodent composition shifted 30% along RDA1 towards the reference as a response to treatment (Figure 5a). Their composition was significantly different between treated and control dunes for semi-fixed dunes. In other words, the ratio between the two rodent species changed significantly as a result of the semi-fixed treatment. Reptiles and beetles shifted along the 1st RDA by 17% and 12% respectively, towards the reference centroid in response to treatment in semi-fixed dunes (Figure 5b,c). The annual plants' centroid shifted away from the reference centroid by 9% along RDA1 (Figure 5d). The centroid shifts detected for these three taxa were not significant in pairwise Permanovas. The 2nd RDA axis for all taxa explained a very low percentage of the variance (0.7-3.3%) compared to the first axis (11.3-49.5%). Distance bars along the 1st RDA axis show the treatment effect (yellow arrow) as a percentage of total distance between control and reference (black dashed); the % is calculated under the scenario that pre-disturbance state was semi-fixed. Significant differences were tested using pairwise Permanova. * and + represent p < 0.05 for treated dunes compared to control and reference dunes, respectively.

Temporal Trends in Composition
To visualize the multivariate changes over time, we plotted the mean scores for each dune type on the first RDA axis (RDA1) from Figures 5-7 as the y-axis, as a function of year along the x-axis ( Figure 8). Rodents on treated semi-fixed dunes were significantly different from their untreated pairs in 2014 (Figure 8a). Beetle composition in treated semi-fixed dunes changed incrementally over time towards the mobile dune centroid, and was eventually significantly different from untreated semifixed dunes in 2016 (Figure 8g). A small shift was seen for reptiles and rodents from untreated to treated semi-fixed dunes in 2016, but it was not significant (Figure 8a,d). These findings could be seen as partial support for hypothesis H4a.
All faunal taxa had significantly different composition in treated semi-fixed dunes and reference (mobile) dunes across all years (Figure 8a,d,g). The same was true for treated fixed dunes ( Figure  8b,e,h). In summary, none of the treated dunes fully reached the reference state in any year for faunal taxa.
For annuals, the semi-fixed treatment had no effect on RDA1 (Figure 8j). Fixed dune treatment appeared to have no effect on RDA1 for any taxa (Figure 8b,e,h,k) supporting the resistance hypothesis H4c . Conversely, disturbed dunes were significantly different from untreated semi-fixed dunes across all years (Figure 8c,f,i) suggesting turnover (H4a). In addition, disturbed dunes were significantly different from mobile dunes in all years for reptiles (Figure 8f), in 2012 and 2015 for rodents ( Figure 8c) and in 2015 for beetles (Figure 8i). In all other years, rodents and beetles were similar to (i.e., not significantly different from) mobile dunes, in RDA1.   Redundancy analysis for (a) rodents, (b) reptiles and (c) beetles in disturbed dunes (off-road vehicle disturbance) compared to semi-fixed untreated (control) and mobile dunes (reference). Distance bars along the 1st RDA axis show the treatment effect (yellow arrow) as a percentage of total distance between control and reference (black dashed); the % is calculated under the scenario that pre-disturbance state was semi-fixed. Significant differences were tested using pairwise Permanova. * and + represent p < 0.05 for treated dunes compared to control and reference dunes, respectively. Rodent composition on control fixed dunes consisted of just one species (G. a. allenbyi) and the second species did not reappear on treated dunes within the timeframe of the study (Figure 6a). Fixed dune treatment produced a minor shift in the intended direction towards the reference for reptiles (2%) and beetles (6%) (Figure 6b,c), while a greater shift (17%) away from mobile dunes was detected for annuals ( Figure 6d). Treated dunes' centroids were not significantly different from controls in any taxa in pairwise Permanovas (see Table S4). The composition of rodents in disturbed dunes shifted from two species found on both mobile and untreated semi-fixed dunes, to one psammophilic species forming an arch-effect in ordination (Figure 7a). Species composition was significantly different to both mobile and un-treated semi-fixed dunes for rodents and reptiles (Figure 7a,b). The centroids for both of these taxa on disturbed dunes were beyond (left of) the centroids of mobile dunes.
The pre-disturbance state for ORV dunes could have been mobile or semi fixed. Under the latter conditions, one could infer a 155% and 261% shift from control to reference along the 1st RDA axis for rodents and reptiles respectively. For beetles, the composition shift along the 1st RDA axis was 88% of the distance between control and reference dunes and was only significantly different from semi-fixed controls in Permanova (Figure 7c). Note that for reptiles, significant differences between disturbed and mobile dunes were found along the primary RDA axis (Figure 7b), while these differences were reflected in the second RDA axis in the overall ordination shown in Figure 4b. This is because the relative difference in composition between disturbed and mobile dunes is much smaller than the difference between mobile and untreated semi-fixed or fixed dune states, irrespective of treatment type.

Temporal Trends in Composition
To visualize the multivariate changes over time, we plotted the mean scores for each dune type on the first RDA axis (RDA1) from Figures 5-7 as the y-axis, as a function of year along the x-axis ( Figure 8). Rodents on treated semi-fixed dunes were significantly different from their untreated pairs in 2014 (Figure 8a). Beetle composition in treated semi-fixed dunes changed incrementally over time towards the mobile dune centroid, and was eventually significantly different from untreated semi-fixed dunes in 2016 (Figure 8g). A small shift was seen for reptiles and rodents from untreated to treated semi-fixed dunes in 2016, but it was not significant (Figure 8a,d). These findings could be seen as partial support for hypothesis H 4 a. Lastly, we conducted a final RDA ordination for each taxon, this time including all dune types across all years ( Figure 9). This allowed us to visualize the relative differences in composition across all dune sites within the study area. Irrespective of treatment, semi-fixed dune composition was more similar to fixed dunes than mobile dunes for rodents and beetles (Figure 9a,c), and relatively similar to mobile dunes compared to fixed dunes, for annuals and reptiles (Figure 9b,d). While disturbed dunes were significantly different to mobile dunes in many samples (Figure 8c,f,i), these dunes were still the most similar to mobile dunes relative to other treatments for all faunal taxa (Figure 9a-c). Figure 10 summarizes the overall findings for abundance, richness and community level responses across all treatment described above. All faunal taxa had significantly different composition in treated semi-fixed dunes and reference (mobile) dunes across all years (Figure 8a,d,g). The same was true for treated fixed dunes (Figure 8b,e,h). In summary, none of the treated dunes fully reached the reference state in any year for faunal taxa.
For annuals, the semi-fixed treatment had no effect on RDA1 (Figure 8j). Fixed dune treatment appeared to have no effect on RDA1 for any taxa (Figure 8b,e,h,k) supporting the resistance hypothesis H 4 c. Conversely, disturbed dunes were significantly different from untreated semi-fixed dunes across all years (Figure 8c,f,i) suggesting turnover (H 4 a). In addition, disturbed dunes were significantly different from mobile dunes in all years for reptiles (Figure 8f), in 2012 and 2015 for rodents ( Figure 8c) and in 2015 for beetles (Figure 8i). In all other years, rodents and beetles were similar to (i.e., not significantly different from) mobile dunes, in RDA1.
Lastly, we conducted a final RDA ordination for each taxon, this time including all dune types across all years (Figure 9). This allowed us to visualize the relative differences in composition across all dune sites within the study area. Irrespective of treatment, semi-fixed dune composition was more similar to fixed dunes than mobile dunes for rodents and beetles (Figure 9a,c), and relatively similar to mobile dunes compared to fixed dunes, for annuals and reptiles (Figure 9b,d). While disturbed dunes were significantly different to mobile dunes in many samples (Figure 8c,f,i), these dunes were still the most similar to mobile dunes relative to other treatments for all faunal taxa (Figure 9a-c).

Figure 10.
A summary of the responses for three dependent parameters to the three independent treatments, across four different taxa. Note some parameters showed no significant change, which could still be the 'desired' response if there was no difference between control and reference dunes.
The comparative strengths of the responses are not depicted.   Figure 10. A summary of the responses for three dependent parameters to the three independent treatments, across four different taxa. Note some parameters showed no significant change, which could still be the 'desired' response if there was no difference between control and reference dunes. The comparative strengths of the responses are not depicted.

Indicator Species
Twenty-eight species were identified as suitable candidates to be indicator species (IS) based on the IndVal analysis, including two rodents, seventeen annuals, nine beetles, and four reptiles ( Table 1). The Indicator Value or IV (a measure of habitat affinity) ranged from 61% to 98% across all taxa. Weighted PPC (wPPC) ranged from 11% for the ground beetle Scarites striatus to 34% for the annual Bromus rigidus. Only one reptile (Stenodactylus sthenodactylus) and five annual species were identified as potential IS for semi-fixed dunes. Most annuals (n = 11) were fixed dune affiliates with high SSI. Only one annual plant species (Cutandia memphitica) was identified as a potential IS for mobile dunes, although it had a relatively low IV of 64% and relatively high wPPC of 18%. The species with the highest and lowest wPPC were selected as IS for further analysis, such that every taxonomic group had one Mobile dune Indicator Species (MIS) and one Fixed dune Indicator Species (FIS) (species in bold in Table 1). The FIS for rodents is actually an indicator for semi-fixed dunes, as identified by the IV analysis, but we refer to it as an FIS for convenience to compare with other taxa since no true FIS was identified. Full model outputs for both MIS and FIS are given in Supplementary material Tables S5 and S6. As expected, the abundances of MIS were much lower (or absent) on untreated semi-fixed dunes compared to the mobile dunes (which is what defines them as good IS with high IV). No change in abundance was visible for the gerbil G. pyramidum (Figure 11a) in treated semi-fixed dunes. A small increase in abundance for the lizard A. scutellatus was also observed in treated semi-fixed dunes, becoming more similar to levels found in mobile dunes, but the increase was not great enough to become significantly different from untreated semi-fixed dunes (Figure 11b). Meanwhile a significant increase in abundance was observed for S. striatus in treated compared to untreated semi-fixed dunes, although it did not reach reference levels (Figure 11c). The annual MIS, C. memphitica had a large range of cover across samples, and cover in treated semi-fixed dunes was not significantly different to mobile or untreated semi-fixed dunes (Figure 11d).  , showing the 25% and 75% percentiles of the measure of abundance (Y axis units e.g., rank abundance/ % cover) on each dune type for rodents, reptiles, beetles and annuals. Dune types were tested for significant differences between treated dunes, their control pair and reference mobile dunes , using mixed effect models. Disturbed dunes were tested against semi-fixed control dune and mobile dunes. Significant results (p < 0.05) are shown for treated dunes compared to *control and + reference dunes (dune type abbreviations refer to the Nizzanim LTER coding system for continuity).
While not significantly different from control, we calculated that both the rodent and reptile MIS abundance increased by 15% of the difference between control and reference levels overall (Table 2i a and b). This can be considered as 15% of the desired response, if the treated dunes were expected to become identical to mobile dunes (H6). For beetles the increase was 27%, and for annuals there was overcompensation in the MIS abundance, with a 114% increase compared to reference ( Table 2). , showing the 25% and 75% percentiles of the measure of abundance (Y axis units e.g., rank abundance/ % cover) on each dune type for rodents, reptiles, beetles and annuals. Dune types were tested for significant differences between treated dunes, their control pair and reference mobile dunes, using mixed effect models. Disturbed dunes were tested against semi-fixed control dune and mobile dunes. Significant results (p < 0.05) are shown for treated dunes compared to *control and + reference dunes (dune type abbreviations refer to the Nizzanim LTER coding system for continuity).
While not significantly different from control, we calculated that both the rodent and reptile MIS abundance increased by 15% of the difference between control and reference levels overall (Table 2i a  and b). This can be considered as 15% of the desired response, if the treated dunes were expected to become identical to mobile dunes (H 6 ). For beetles the increase was 27%, and for annuals there was overcompensation in the MIS abundance, with a 114% increase compared to reference (Table 2). Table 2. Effectiveness of treated dunes on abundance of Mobile and Fixed-dune Indicator Species (IS). Calculated as the change (either increase or decrease) in abundance (faunal taxa) or % cover (annuals) from control to treated dune, as a percentage of the difference in abundance between control and reference dunes for each species.

Indicators
Taxa The abundances of the MIS G. pyramidum and A. scutellatus on disturbed dunes were also significantly higher than their abundances in semi-fixed control dunes and were similar to the abundances found on mobile dunes (Figure 11a,b), increasing by 126% and 68% of the expected abundance respectively ( Table 2). The abundance of S. striatus was significantly higher on disturbed dunes compared to semi-fixed control dunes (42% increase), but it was still significantly lower than on mobile dunes (Figure 11c). None of the MIS responded to fixed dune treatments rejecting H 6 ; they remained absent on treated fixed dunes (Figure 11a-d).
FIS was not found to show any significant response to fixed dune treatment in any taxa (Figure 11e-h) leading us to reject H 6 . Although the abundance of G. a. allenbyi significantly decreased in treated semi-fixed dunes compared to their untreated pair, the response was not enough to become similar to mobile dunes (Figure 11e). Lastly, disturbed dunes had a significantly lower abundance of this rodent than both mobile and untreated semi-fixed dunes, while the beetle FIS was only found to be significantly lower than untreated semi-fixed dunes, being almost entirely absent from mobile dunes (Figure 11g).
Mean abundance of rodent and beetle FIS declined by 30% and 33% of the desired reduction respectively ( Table 2). Annual FIS % cover declined by 21%, and there was no change in fixed dunes for reptile FIS. In addition, mean abundance decreased for all FIS in semi-fixed treated dunes by 30%, 36%, 30% and 42% in rodents, reptiles, beetles and annuals respectively, towards the reference.
We plotted the abundance of each species (mean for each dune type) over time to show populations fluctuations ( Figure 12). These abundances were tested yearly in a pairwise permuted distribution test. The significance of these tests should be interpreted with caution due to the small sample sizes and multiple testing. Nevertheless, some interesting trends were apparent across the different taxa.
All MIS abundances fluctuated stochastically on the mobile dunes over time (Figure 12a-d). They were all observed with relatively low to zero abundance on semi-fixed control dunes, with the exception of the lizard A. scutellatus, with intermediate abundance in semi-fixed dunes. All MIS were effectively absent from fixed control dunes. In semi-fixed treated dunes, the abundances of both the rodent G. pyramidum and the beetle S. striatus were similar to untreated semi-fixed dunes (i.e., low abundance) and were significantly different from AC, until 2016. In the final year of monitoring, the abundances of both these species increased and they became similar (not significantly different from) mobile dunes and significantly different from untreated semi-fixed dunes (Figure 12a,c). The same was partially true of the annual C. memphitica, which was similar to untreated semi-fixed dunes until the final year 2016, but was not consistent across years in comparison to mobile dunes (Figure 12d). The lizard A. scutellatus showed a small but similar trend in the final year of monitoring in treated semi-fixed dunes, but it remained significantly different from mobile dunes for all years, and similar to untreated semi-fixed dunes in all years except 2013 (Figure 12b). Indicator Species (FIS) on different dune types for rodents, reptiles, beetles, annuals (see icons by row). Dune types were tested for significant differences between treated dunes, their un-treated control pair and the reference mobile dunes each year. * and + represent p<0.05 for treated dunes compared to control and reference dunes, respectively. Not all significant results are shown given for intelligibility. See Supplementary Tables S7 and S8 for full results.
In the disturbed dunes, all FIS were effectively absent, with the exception of G. a. allenbyi, whose abundance was low to medium on both disturbed and mobile dunes in most years (Figure 12e). This species was not detected in disturbed dunes in 2013, while it was still observed in mobile dunes. The abundance of this species in untreated semi-fixed dunes was significantly higher than in disturbed dunes for all years. A. schreiberi was present occasionally in low abundances on mobile and untreated semi-fixed dunes over the course of the study, but remained absent in disturbed dunes in all years monitored (Figure 12f). G. serrator was effectively absent on both disturbed and mobile dunes in all years, and only occasionally appeared in low abundances on untreated semi-fixed dunes. Thus, abundance of this FIS in disturbed dunes was found to be similar to both mobile and untreated semifixed dunes in all years (Figure 12g).

Discussion:
Year G e rb illu s p y ra m id u m G e rb illu s a n d e rso n i a l l e n byi Dune types were tested for significant differences between treated dunes, their un-treated control pair and the reference mobile dunes each year. * and + represent p<0.05 for treated dunes compared to control and reference dunes, respectively. Not all significant results are shown given for intelligibility. See Supplementary Tables S7 and S8 for full results.
In the fixed dunes, none of the MIS were found to have differences in abundance between treated and un-treated fixed dunes (i.e., they remained absent or in very low abundance), and all the faunal MIS remained significantly different from mobile dunes in all years (Figure 12a-c).
In the disturbed dunes, G. pyramidum had similar abundances to mobile dunes in all years except 2014 (Figure 12a). The abundance of S. striatus was significantly lower in disturbed than in mobile dunes in all years except 2014, although there did appear to be a trend of increasing abundance in disturbed dunes over time (Figure 12c). In addition, the latter species was observed in disturbed dunes in all monitored years, while it remained absent from un-treated semi-fixed dunes throughout the study. In fact, both G. pyramidum and S. striatus had significantly higher abundances in disturbed dunes, compared to untreated semi-fixed dunes in all years (Table S7). Abundance of A. scutellatus on disturbed dunes was most similar to mobile dunes in all years except 2015, and was significantly different to untreated semi-fixed dunes in 2013 and 2014 (Figure 12b). Next, we considered the effect of each treatment on FIS (Figure 12e-h). In semi-fixed dunes, the rodent G. a. allenbyi initially had similar abundances in treated and untreated semi-fixed dunes, which were significantly higher than in mobile dunes, but later in the study, the abundance of this species declined and became more similar to mobile dunes and less similar to untreated semi-fixed dunes (Figure 12e). The abundance of the lizard A. schreiberi was similar across mobile, treated and untreated semi-fixed dunes (i.e., low to medium abundance), except in 2012 where it was absent in mobile dunes, but remained present in both treated and untreated semi-fixed dunes (Figure 12f). The abundance of ground beetle G. serrator was initially similar (low abundance) in mobile, treated and untreated semi-fixed dunes. Then in later years, abundances in treated semi-fixed dunes became different to mobile dunes, but remained similar to untreated semi-fixed dunes (both increasing slightly) in 2015 and 2016 (Figure 12g). The annual plant B. rigidus had generally very low cover in all dune types (Figure 12h). The cover in treated semi-fixed dunes was similar to mobile dunes in all years (i.e., absent) and similar to untreated semi-fixed dunes in all years except 2014, when some presence was observed on untreated semi-fixed dunes, while remaining absent from treated semi-fixed dunes.
In the treated fixed dunes, abundance of G. a. allenbyi was initially different to mobile dunes but became similar in the later years of the study (Figure 12e). However, abundance of this species in DC also showed declines until 2009, and abundance in treated and untreated fixed dunes was similar in all years except 2005 and 2008. Unfortunately, data is lacking for the period between 2009-2016 on treated fixed dunes, and trends in the interim years are unknown. The other FIS, A. schreiberi, G. serrator, and B. rigidus had similar abundances on treated and untreated fixed dunes throughout all years, despite stochasticity (Figure 12f-h). A marked increase in cover was observed in the last year of the study for B. rigidus, but this was detected in both treated and untreated fixed dunes, so overall, no difference was found. Both G. serrator and B. rigidus were absent from mobile dunes in all years. All three FIS were significantly different from mobile dunes in all years (Figure 12e-g), excluding years where B. rigidus was absent in both treated and untreated fixed dunes and mobile dunes (Figure 12h).
In the disturbed dunes, all FIS were effectively absent, with the exception of G. a. allenbyi, whose abundance was low to medium on both disturbed and mobile dunes in most years (Figure 12e). This species was not detected in disturbed dunes in 2013, while it was still observed in mobile dunes. The abundance of this species in untreated semi-fixed dunes was significantly higher than in disturbed dunes for all years. A. schreiberi was present occasionally in low abundances on mobile and untreated semi-fixed dunes over the course of the study, but remained absent in disturbed dunes in all years monitored (Figure 12f). G. serrator was effectively absent on both disturbed and mobile dunes in all years, and only occasionally appeared in low abundances on untreated semi-fixed dunes. Thus, abundance of this FIS in disturbed dunes was found to be similar to both mobile and untreated semi-fixed dunes in all years (Figure 12g).

Discussion
Coastal dune habitats have declined globally over the last several decades due to urbanization. Within remaining dune systems, substantial dune fixation has resulted in further losses of mobile dunes with negative impacts on their associated species [15,19,25,27,36,49,78,94]. Some studies suggest vegetation removal can initially promote habitat heterogeneity, and increase availability of suitable habitats for psammophile, xeric and endemic mobile dune species, but longer term responses are generally unknown [15,24,28,40,54]. We investigated the temporal trends of four taxonomic groups (including fauna) to determine the effect of vegetation removal on dune assemblages over a 12-year period. Detecting community level responses to restoration can be difficult. Indeed, two studies in Californian coastal dunes were unable to show significant assemblage level responses for terrestrial arthropods [32] or plants [24] despite intensive vegetation removal. Our study found only marginal restoration effects in most cases.
This large-scale project across multiple taxa and multiple treatments meant that our samples sizes were relatively small, and some caution is needed to interpret the results. Nevertheless, some inferred conclusions can be made. Comparison of control dunes only supported H 1 for annual plants (% cover increases with perennial cover) and H 2 for annuals and reptiles (richness increases with perennial cover). Composition level parameters showed only small shifts away from the control dune, suggesting high levels of resistance at the community level. Although many of these changes were not statistically significant, they were usually in the direction desired by the manipulation, towards mobile dune status. Semi-fixed dune responses were more pronounced than fixed dune responses. The ORV-disturbed dunes shared more aspects with mobile dunes for faunal taxa compared to both semi-fixed and fixed treated dunes, while annuals were effectively absent. Most mobile dune indicator species also showed some trends in the intended direction in response to most treatment types suggesting that treatment did increase habitat availability for these species, while not removing suitable habitat for fixed dune species. Our results suggest that more pronounced responses might occur given a longer time frame, and that press disturbance is more effective than pulse disturbance for restoring mobile dune habitat.
Overall, abundance and richness measures were not appropriate indicators for treatment success (since there were few differences between control dune states). Richness did correlate with the gradients of increasing PPC across control dune dunes, but the effects of treatment detected for this parameter were marginal. Both abundance and richness are important to monitor in restoration projects, but they are more appropriate for assessing whether any unexpected artefacts resulted from the intervention. Composition is known to a better parameter for detecting changes [26,28]. The indicator species that were selected also were highly suitable for detecting changes. Both composition and indicator species showed a temporal trend in treated semi-fixed dunes, away from the control dunes towards the reference dunes over time; this trend was more consistent in the faunal taxa. Over a period of 12 years, we found that the response at the assemblage level was between −9% to 30% in semi-fixed treatment, compared to just −17% to 6 % in fixed dunes.

Effect of Treatment on Semi-Fixed Dunes
Establishing dune dynamics on the wind facing slope can lead to a rejuvenation of the landscape through reactivation of sand, leading to burial of any newly developing vegetation [31]. This intervention enables pioneer species found in mobile dunes to recolonize [4]. Rodents responded well to treatment in semi-fixed dunes with a significant shift in composition detected, achieving 30% of the desired Euclidean ordination distance towards reference, as expected by Hypothesis H 4 . In addition, though the change was not significant for reptiles or beetles, the response to treatment was consistently in the desired direction towards mobile-dune composition. For annuals the shift was minor, but away from the desired composition, becoming less similar to mobile dunes, and H 4 was rejected.
At the temporal scale, rodents showed significant differences in some years. The change in reptile and beetle composition was more pronounced in the final years of the study, and became significant for beetles in the last year. Similar temporal responses were also detected at the indicator species (IS) level, particularly for "mobile-dune indicator species" (MIS). These results indicate some support for Hypothesis H 4 a, where the composition turnover is due to recolonisation by MIS.
The IS were selected objectively based first on their affinity to fixed and mobile dunes, and secondly as the species with the lowest and highest weighted PPC. Three MIS (all except the reptile), as well as the rodent "fixed-dune indicator species" (FIS) appeared to respond desirably in the direction intended for restoration. Given that the beetle S. striatus was almost entirely absent from control dunes and increased over time in treated dunes, we can infer that its increased abundance was due to colonization rather than competitive release. The small increase detected for rodents and reptiles was not significant when looking at averages across years, but the gerbil G. pyramidum did respond favorably over time, showing significant increases in the final year of the study. This increase may be due to the overall decline in abundance of the subordinate G. a. allenbyi populations, since these two species compete sympatrically at the patch scale [95][96][97]. Further investigation is needed, since the trends of the two species were not consistent over time. The annual MIS C. memphitica also had relatively high cover in some years on treated dunes, while remaining much less prevalent in control dunes. This was most likely a response to competitive release following the removal of perennial vegetation. Over time, this species was erratic in both mobile and treated dunes and the temporal trends were not synchronized, making it a poor IS, as expected due to the relatively low IV.
Given similarities between control and reference dunes in terms of reptile and beetle species richness, there was no expectation for treated dunes to elicit any response. However, control and reference dunes were different for annual plants, so for this taxon, a lack of response can be seen as lack of desired treatment effect. Only rodents increased in richness as expected under Hypothesis H 3 b (from one to two species) and H 3 a was rejected for all taxa in semi-fixed treated dunes. These results highlight the importance of examining composition rather than univariate measures when studying restoration effects on biodiversity.
The composition shift in semi-fixed treated dunes may continue towards a mobile state with more species responding positively, but without a press disturbance, we may observe a return to the control state. In the Netherlands, aeolian activity increased after an intervention, up to a climax in the third year, and then started to decline because of vegetation development on bare spots [4]. Arens et al. [40] suggested that follow-up removal of roots is needed for three to five years and it is not clear whether the dunes will remain mobile in the long term.

Effect of Fixed-Dune Treatment
A recent review of coastal restoration projects revealed that fixed dunes are the most frequently restored dune type (54%), followed by semi-mobile dunes (30%), mobile dunes (27%) and dune slacks (16%; [2]. However, the most frequently used restoration mechanism was revegetation (42%), while destabilization or vegetation removal only accounted for less than 5% of restoration projects. Removal methods varied among removal experiments.
The extensive and intrusive removal of perennial vegetation in a grid formation on fixed dunes had virtually no effect on most of the parameters monitored, rejecting hypothesis H 3 b, H 4 a, and H 6 for this treatment type. The only consistent effect was a reduction in the average abundance of all taxa (H 3 a); richness and composition were not affected. That said, at the temporal scale, some weak responses of reptiles were observed both for the IS and for composition in the first two years after treatment, suggesting some support for resilience under Hypothesis H 4 b. All other taxa showed very limited response to fixed dune treatment, indicating strong resistance to change (supporting Hypothesis H 4 c). Detachment from mobile dune assemblages may explain this lack of response. In addition, the dunes from which vegetation was removed are surrounded with dense interdune vegetation, including tall sycamore trees, which reduce wind speed, preventing sand migration and potential recolonization from neighboring areas. Thus, sand deposition in the treated dunes was minimal, and sand enriched with nutrients and organic matter remained present after the vegetation removal. Throughout Europe, several coastal dune management schemes involving vegetation removal achieved partial success in restoring dune mobility, from a geomorphological perspective [6,15,31,40].
Evidence suggests that soil processes in fixed dunes have changed the soil irreversibly [19]. In experimental dune restorations, removal of vegetation and topsoil in the slacks and inner dunes often fail to restore aeolian activity due to the groundwater table and the usually moist conditions at the surface, which prevent saltation [4]. In Mediterranean coasts, wind-related variables may only be of minor importance compared to soil properties [98]. Encroachment and intended establishment by grasses, shrubs and woody plants have accelerated in response to high rates of nitrogen deposition [33,[99][100][101][102]. Top-soil inversion has become a popular method for restoration of inland dunes in recent years, where topsoil is buried under a layer of subsoil, so that the original layers are intact but their position in the profile is changed [103,104]. Results indicate that inversion can successfully lower surface soil fertility and reduce competition, creating a nutrient-poor space for the establishment of early-successional stage species in calcareous grasslands and sand dunes [104,105]. Ploughing and topsoil inversion enhanced lichens and annual plants, lowered organic matter and increased rabbit activity compared to control plots in grass heath habitats [106,107]. Thus, topsoil removal within the grid formation applied on Nizzanim fixed dunes may have created a better response at the community level.

Effect of ORV Disturbance
Off-road vehicle (ORV) disturbance is an anthropogenic press disturbance that is increasingly threatening coastal ecosystems [108,109]. In general, recreational disturbances result in a decrease of species diversity [19,[110][111][112]. Motorbikes, 4 × 4 vehicles and quad bikes are driven on beaches and dunes for recreation purposes [113]. These vehicles can cause more widespread damage than human trampling [109,114,115], another disturbance that is a widely reported issue for fragile coastal biodiversity [110][111][112][115][116][117]. The negative impacts ORVs cause include damage of the physical properties and stability of the substrate, destruction of vegetation, and disturbing, injuring or killing fauna [108,114,[118][119][120][121]. In California, dunes exposed to ORV disturbance showed significant decline in Coleoptera diversity and species evenness [122].
The effect of disturbance on species composition depends on whether the disturbance exacerbates or reduces environmental harshness, and the conditions that favor specialization (Attum 2006). Since recreational driving in the dunes reduces perennial shrub cover and facilitates sand remobilization, this disturbance could theoretically benefit dunes. The situation is complex because we were not able to quantify the start date, frequency, or intensity of vehicular usage. It is speculated that dunes used by ORVs are most likely to have been mobile dunes, but we do not know what trends would have occurred naturally without disturbance over this time period. However, by comparing the disturbed dunes with both mobile and control (un-treated) semi-fixed dunes, we were able to quantify the outcomes of ORV disturbance on species composition. While disturbed dunes were significantly different to mobile dunes in many samples, these dunes were still the most similar to mobile dunes relative to other treatments for all faunal taxa. Nevertheless, there was undoubtedly a negative impact of ORV disturbance in Nizzanim reserve; firstly, annuals were almost entirely absent due to the intensity of disturbance preventing seedling establishment. Secondly, a significant reduction in total abundance and overall richness of reptile species was seen compared to both mobile and semi-fixed states.
Reptiles are known to be particularly sensitive to the presence of vehicular disturbance, affecting the relative abundance and proportional use of preferred microhabitats in coastal systems of Argentina [123]. Rocha and Bergallo [124] found that a gradual reduction in beach vegetation over a decade was accompanied by population declines of a lizard in coastal Brazil. Luckenbach & Bury [120] reported heavy impacts of ORVs on lizards, mammals and arthropods in California. We also found that the populations of the ground beetle S. striatus were significantly lower than on mobile dunes.
If the state of disturbed dunes prior to disturbance was more similar to semi-fixed dunes, then the ORV disturbance achieved a much greater shift in composition on these dunes (towards reference) than either of the pulse treatments achieved. Indeed, if the restoration goal was to increase the abundance of MIS, the response to ORV disturbance in the study site could be perceived as positive.
The abundance of the desert gerbil G. pyramidum was significantly higher in ORV disturbed dunes, compared with both semi-fixed control and mobile dunes. Thus, the overall rodent assemblage was almost entirely composed of this species. Plausibly, the 'hypermobile' conditions created are only inhabitable by highly adapted sand specialist species. Both gerbils are found in sandy habitats, but G. pyramidum is better adapted to exposed, harsh conditions, while G. a. allenbyi needs some shrub cover for protection [97,125]. Population levels of the beetle S. striatus were also significantly higher in ORV dunes than in semi-fixed dunes, although still somewhat lower than numbers found in mobile dunes. This species is also highly adapted to sandy desert conditions. Compared to untreated semi-fixed dunes, ORV usage resulted in an increase in highly adapted MIS, but overall alpha diversity was reduced.
In a review of coastal restoration projects [19] only two studies did not find a significant impact of recreation disturbance, and both were conducted in Israel's Mediterranean dunes in other coastal reserves [114,115]. Despite the differences from mobile dunes, the overall composition in ORV-disturbed dunes in Nizzanim remained the most similar to the reference mobile dunes, compared with the two press disturbances for all three faunal taxa. Thus, we found evidence supporting Hypothesis H 5 that press disturbance is more effective than a single pulse. That said, richness in these dunes was generally lower. Therefore, widespread, pervasive ORV disturbance in NDNR will likely result in an overall decrease in γ-diversity. The LTER experimental design was not balanced across dune types or years, therefore these comparisons should be taken as inferred rather than proven.

Implications for Conservation Management
Disturbance is a fundamental part of natural coastal dune processes, and encroachment by grasses, shrubs and woody plants has most likely occurred due to the removal of this process [4,33]. Geomorphologists increasingly recommend that remobilization of inner dunes should consider the dune-beach system as a whole, arguing for an integrated dynamic approach that recreates the natural disturbances such as wind erosion and sand burial [15,24,40,126]. Our results provide support to the growing body of evidence that restoration attempts should promote, not diminish, the degree of disturbance in coastal dune systems. e.g., [31,38,49].
When comparing the pulse treatments, vegetation removal of the wind-facing slopes of semi-fixed dunes appeared to be more effective than grid formation removal of fixed dunes, but since the experiments are not equivalent this can only be inferred. An initial pulse disturbance is likely needed to remove perennial vegetation, but continued press disturbance appeared to more closely match conditions found on mobile dunes, at least for faunal taxa. While the unregulated use of ORVs as a method for maintaining dune mobility is certainly not recommended, previous evidence has demonstrated the benefit of press disturbances to plant and arthropod species [15,31,49]. ORV disturbance has severe impacts on annual plant species, so other forms of press disturbance such as grazing could be more beneficial.
Grazing has been successful in increasing both plant and arthropod diversity on coastal dunes [49,127]. However, domestic grazers may affect the composition (of vegetation) differently than wild grazers [128]. A review of the effects of grazing by domestic livestock showed both positive or negative impact on diversity, often increasing diversity at the landscape level but reducing alpha diversity at the very local scale [19].
Depending on the reference species, grazing can have a positive or negative impact on richness and composition [34,127,129,130]. Ground beetles and reptiles are often negatively affected by management practices such as grazing and trampling [129,131,132]. That said, populations of A. scutellatus or S. striatus on disturbed dunes were not as low as in semi-fixed control dunes, suggesting that even recreational driving is less damaging than passive management (allowing shrub encroachment) for these species. Attum [133] showed that desert reptile specialists in the Egyptian desert did not respond favorably to heavy grazing in inland sand dunes, despite predictions that they may fare better than generalist species. Another study found arthropod assemblages were more diverse in response to disturbance by trampling, blowouts and burning in Danish dune grasslands [44]. Grazing management applied in coastal dune grasslands of France and Belgium was also able to control highly competitive plant species, and may be beneficial in maintaining mobile assemblages [134]. Grazing could be introduced either through re-wilding of the ecosystem or in the form of managed domestic grazers [49,135].
Mountain gazelles (Gazella gazelle gazella), the only mammalian grazers in NDNR, are unlikely to provide an adequate grazing impact in the near future [136]. Attempts to introduce camel grazing to Nizzanim on semi-fixed dunes had little effect on the perennial vegetation [137]. Other forms of grazing or press disturbances should be investigated. Reintroduction of small livestock herds (sheep and goats) are currently being trialed with Bedouin shepherds, and should be monitored carefully before further recommendations can be made with regarding to grazing.

Conclusions
Conservation managers must carefully consider their priorities when implementing coastal restoration actions. Is the conservation goal a specific species that is regionally or locally endangered, or is the goal to preserve mobile dune communities, or increase alpha or gamma diversity? These decisions affect the desired actions and restoration strategy. The scale and extent of the introduced disturbances for management purposes must be carefully planned to prevent erosion and detrimental effects and should be tailored the specific goals [2].
The methods used to identify IS (IV and wPPC) appeared to be a good choice for selecting species, and several of these species appeared to respond favorably to treatments, in particular the MIS. The MIS are of particular interest for restoration references, since they are by definition found only on mobile dunes and are therefore most likely to decline across the reserve as a result of widespread dune fixation. Therefore, the positive albeit weak response of these species is encouraging. This research highlights the potential use of IS in assessing the impact of restoration projects, and emphasizes the importance of long-term studies in restoration projects.
Having longer time-series for pre-restored sites would be beneficial [138]. We strongly encourage LTER sites monitoring biodiversity trends to adopt a regularized temporal sampling regime that enables powerful and simpler statistical analyses [73,139]. That said, empirical ecological research projects in natural conditions are often unable to achieve perfect experimental design due to historical and external factors affecting accessibility. Balanced sampling across multiple taxa in a long-term study is invariably difficult for empirical data [140]. The approach to concurrently monitor treated, control and reference dunes allowed us to make effective interpretations, despite the unbalanced design. We hope the statistical approaches we adopted for unbalanced data will be useful for other researchers and practitioners.
Shrub removal experiments allowed us to investigate the structural component of dune states rather than the biological impact of the primary productivity. It is apparent that the shrubs themselves determine biodiversity trends rather than just providing biotic resources. Each taxon responded differently to removal, and each type of treatment had different effects. Restoration manipulations were only partially successful in Nizzanim LTER, but rodents seemed to respond most favorably. With respect to ecosystem functioning, "long-term" can mean anything from decades to centuries, so realistic expectations about restoration of sand dunes could take many years [141]. Over a longer timeframe, the response at the composition level may become more pronounced. Overall, it was apparent that in the absence of any disturbance, gamma diversity of coastal dunes communities would be lost as dune fixation increases.
It is likely that more 'natural' disturbances such as controlled grazing could be more effective and less damaging than mechanical ORV disturbances. Indicator species can represent assemblage level responses, although their response was often more pronounced than assemblages. Finally, it should be remembered that long-term assessment and monitoring of restored coastal dunes should include multiple taxa, and should examine changes in composition rather than abundance and richness alone.
Supplementary Materials: The following are available online at http://www.mdpi.com/2076-3417/10/7/2310/s1, Table S1: Experimental design for the long-term ecological research project in Nizzanim, showing the sampling time (years) for each taxon, Table S2: Poisson (GLMM_PQL) and Linear (LME) Mixed effect Models of community abundance ( Figure 2) for each taxon, Table S3: Poisson (GLMM_PQL) and Linear (LME) Mixed effect Models of species richness (Figure 3) for each taxon, Table S4: Pairwise Permanova of species composition among dune states. Each treatment type was tested separately, with some data overlap (e.g., mobile and Untreated Semi-Fixed dunes are used for both treated Semi-Fixed and Disturbed dune comparisons). Bold rows are of interest in terms of meeting the expectations for responses to treatment, Table S5: Poisson (GLMM_PQL) and Linear (LME) Mixed effect Models of abundance for each Mobile-dune Indicator Species (MIS). Results were used for Figure 11a