Long transient response of vegetation dynamics after four millennia of anthropogenic impacts in an island ecosystem

Abstract Agents of global change commonly have a higher impact on island ecosystem dynamics. In the Mediterranean region, those dynamics have historically been influenced by anthropogenic impacts, for example, the introduction of invasive species and overharvesting of resources. Here, we analysed the spatio‐temporal dynamics of vegetation in sa Dragonera island, which experienced a large environmental change ca. 4000 years ago by the arrival of humans. Anthropogenic impacts, such as herbivory by goats and over‐logging, ended in the 1970s, while in 2011 the site became the largest Mediterranean island where rats were eradicated. Invasive rats and goats played the ecological role of two endemic species, the cave goat and the giant dormouse, which inhabited the island for more than 5 million years and were rapidly extinct by humans. We used Landsat imagery to explore NDVI as a proxy of vegetation productivity over the years 1984–2021, orthophotos to assess changes in land and vegetation covers and historical plant inventories to study the dynamics in plant diversity. Results showed that those indicators steadily increased both in spring and in summer, while the noise around the trends was partially explained by climate variability. The regime shifts in the temporal dynamics of vegetation productivity suggested a transient from a perturbed to a non‐perturbed stable state. Trends in successional dynamics, spatial self‐organization and plant diversity also showed the same type of transient dynamics. Historical perturbations related to harvesting (mainly the synergies between goat browsing, burning and forest over‐logging) were more important than rat eradication or the influence of climate to explain the vegetation dynamics. Our study shows the transient nature of this small island ecosystem after 4000 years of perturbations and its current path towards vegetation dynamics more controlled by ecological interactions lacking large herbivores and omnivores, drought dynamics and the carrying capacity of the island.


| INTRODUC TI ON
The influence of anthropogenic activities on ecosystems has been particularly intense in the Mediterranean region at least since the second millennium BC (Blondel et al., 2010;Médail, 2017). A temperate and relatively stable climate allowed high human population density, which brought about long-term harvesting of natural resources (e.g. wood and wildlife), habitat loss for cultivation and rapid processes of domestication of plant and animal species, among other effects. Furthermore, a closed sea with calm conditions and relatively close seashores allowed humans to sail for warfare and trade interests, which also moved organisms (e.g. plants, animals and fungi) intentionally or unintentionally throughout the entire region (Blondel, 2006;Coll et al., 2010). This intense maritime traffic particularly affected the biota of Mediterranean islands and its ecological processes, which evolved mainly by climate variability and ecological interactions since the end of the Mesiniense crisis in the upper Miocene, ca. 5.3 million years ago (Alcover et al., 2000;Mas et al., 2018;Masini et al., 2008). As recorded at the global level, islands have been particularly affected by the arrival of humans, due to their fragility to harvesting in spatially reduced and simpler ecosystems, and the biogeographical barrier that sea represents (i.e. slower colonization processes and the hampering of demographic rescue from source populations) (Blackburn et al., 2004;Wood et al., 2017).
To fully understand the functioning of island ecosystems and their processes as we see them nowadays, it is necessary to explore the composition and dynamics of ecological communities over deep time (Carrión et al., 2010;Fraser et al., 2021). Palaeontological records from the Balearic archipelago, in the western Mediterranean, show that settlement of humans ca. 4000 years BP represented a transition to a new ecological state, in which species and communities lacking evolutionary responses to human harvesting and the arrival of new species disappeared or kept relict populations   (Figure 1a). Most endemic fauna of the Balearic islands was extinct and many surviving species disappeared from the large islands and remained refugees in the satellite islets after the arrival of humans (Alcover, 2010;Bover & Alcover, 2003;Mas et al., 2018). The Balearic fauna since the onset of the Pliocene 5.3 million years ago lacked large terrestrial predators, such as carnivores, but some islands (e.g. Mallorca) held two species relevant for the present study: the dwarf herbivore Balearic Islands cave goat (Myotragus balearicus) and the omnivore rodent Majorcan giant dormouse (Hypnomys morpheus) (Alcover et al., 2000). These endemic species likely had enlarged ecological niches due to the lack of intraguild competitors as they were the only ungulate and rodent of the Majorcan fauna, and their populations were top-down regulated by aerial predators known to occur here at that time, that is, golden eagles (Aquila chrysaetos) and the giant barn owl (Tyto balearica) respectively ( Figure 1a). Bottom-up regulation in this simplified food-web was completed by shrublands, especially large patches of Balearic boxwood (Buxus balearica) browsed by the cave goat, which was adapted to the poisonous toxin of the shrub (Alcover et al., 2000;Alcover, Perez-Obiol, et al., 1999). Palaeodietary studies on cave goats and giant dormice showed that, together with climate conditions, they influenced vegetation dynamics and the structure of plant communities (Alcover, Perez-Obiol, et al., 1999;Burjachs et al., 2017;Hautier et al., 2009).
For understanding the changes in plant communities we may observe in recent decades in the small island of sa Dragonera, just 900 m off Mallorca, it is essential to explore the dynamics before and after the arrival of humans. Majorcan palynological records show a sudden shift in vegetation composition caused by a warming period around ca. 4900 years BP, where a community dominated by junipers Juniperus, deciduous trees (mainly Quercus spp.), hazel trees (Corylus spp.) and Balearic boxwood changed towards a maquiatype landscape with a higher prevalence of wild olive (Olea), pines (Pinus) and evergreen oaks (Quercus spp.) (Burjachs et al., 2017). Ca. three millennia later, humans moved here, among others, two species that were particularly successful to establish large populations: domestic goats (Capra hircus) and brown rats (Rattus rattus), which were ecologically close to the original endemic cave goat and giant dormouse respectively. After the introduction of goats and sheep, humans likely burnt forest patches, including Balearic boxwood, since domestic herbivores have not developed a defence against its toxin (Alcover et al., 2000). Humans also logged timber for wood and fuel and for opening space for cultivation. Harvesting of forests and grazing by free-ranging domestic herbivores were synergetic with the abovementioned influence of climate change and determined a maquia-type landscape, which shaped a new ecological state that has since then framed the dynamics of plant communities until recently. A new transition to a different state of vegetation dynamics is expected since the cessation of intense harvesting of natural resources, especially in small islands of the Balearic archipelago around the mid-1970s (Carrión et al., 2010;Picornell-Gelabert et al., 2021).
Rural abandonment and the arrival of conservation laws a few years later allowed managers to protect these small ecosystems from new anthropogenic impacts, mainly habitat loss by touristic housing development. Sa Dragonera Island, with a surface of 276 ha, is an example of such an ecosystem: plant communities have been shaped over centuries by intense harvesting practices, mainly grazing by goats, omnivore foraging by rats, shrub harvesting and forest logging ( Figure 1). Furthermore, in 2011 sa Dragonera became the largest Mediterranean island where brown rats, mice and rabbits were eradicated in a single-day operation (Ibañez-Álvarez et al., 2022). Rat occupation has been a global ecological driver in the whole region: 99% of Mediterranean islands in the western basin with a >30 ha surface are inhabited by rats, which have been held responsible for severe impacts on native biota Latorre et al., 2013;Martin et al., 2000;Ruffino & Vidal, 2010;Traveset et al., 2009).
Sa Dragonera represents a suitable ecosystem to disentangle the role of climate and anthropogenic impacts on plant communities in the long term, and to study the transient nature of a return path towards a state without those impacts ( Figure 1b). Specifically, we studied several features of vegetation dynamics at a landscape and community scales. We first explored the temporal dynamics of vegetation productivity, the influence of climate and the occurrence of regime shifts caused by rat eradication using the NDVI index (normalized difference vegetation index) from Landsat satellite images over the last four decades . NDVI is a simple indicator and powerful tool used for assessing spatio-temporal changes in vegetation distribution, productivity and dynamics (Pettorelli, 2013).
Second, we studied the dynamics of land cover , which has proven to be suitable for assessing decadal-scale structural changes in vegetation caused by climate and anthropogenic impacts (Cousins, 2001;Silapaswan et al., 2001). Third, we assessed the spatio-temporal dynamics in plant cover and spatial pattern focusing on changes caused by rat eradication and considering the degree of self-organization, which is a pattern typical of non-perturbed ecosystems (Berdugo et al., 2017;Kéfi et al., 2007Kéfi et al., , 2011. Finally, we used all plant inventories collected in the field during 1921-2021 to study the dynamics in plant species richness. Plant diversity may show long transients lasting many decades after perturbations due to the dynamics of ecosystem interactions (Fournier et al., 2020).
Most studies dealing with vegetation dynamics over time analyse a single feature of vegetation dynamics (e.g. productivity) using a specific methodology. In our study, we applied the most up-to-date reliable methods for exploring the dynamics of vegetation for each of its features (Ahmed et al., 2013;Berdugo et al., 2017;Ibañez-Álvarez et al., 2022;Zhang et al., 2018).
Despite most of the available data started ca. 10 years after the cessation of the abovementioned anthropogenic perturbations (Figure 1), we expected that those data would still show an increase in productivity, cover and self-organization, that is, the spatial pattern through the biotic interactions of facilitation and competition that operate in the successional process. We also expected an increase in forest land cover to the detriment of shrublands, and an increase in plant diversity over the last century. Furthermore, all study F I G U R E 1 (a) Schematic reconstruction of the ecosystem state in sa Dragonera since the establishment of the endemic fauna over the Pliocene, ca. 5.3 million years ago and anthropogenic impacts since the arrival of humans, ca. 4 millennia ago. Most palaeontological records come from the island of Mallorca, but due to its close proximity, we assumed that the endemic fauna was similar in sa Dragonera. Extinctions occurring in the small sa Dragonera were likely followed by colonizations from Mallorca, especially during glaciation events and flotsam events. Red arrows show extinctions, grey arrows refer to citations of impacts or their ending and black arrows show the four crucial events to interpret NDVI vegetation dynamics during 1984-2021 (blue arrow). Sparse old data on vegetation features were available for 1921 and 1956 (blue dots). Years correspond to approximate or specific dates recorded in archaeological data and historical and contemporary chronicles (more details are provided in Supplementary Material) (b) Schematic reconstruction of vegetation dynamics over the same temporal window representing the three stable states: herbivory regulated by predation before the arrival of humans (K 1 ), harvesting and not regulated herbivory (K 2 ), and without herbivory and harvesting (K 3 ). Yellow lines represent the expected NDVI major transient dynamics.
vegetation parameters (productivity, cover, spatial self-organization and diversity) should have increased at a higher rate after rat eradication due to the release of detrimental effects that are commonly attributed to rats in island ecosystems (Howald et al., 2007;Traveset et al., 2009). Finally, we tested the hypothesis that transient dynamics have occurred for vegetation dynamics (Figure 1b). Transient would be consistent with a return towards a new ecosystem structure after cessation of perturbations (Hastings et al., 2018).

| Study sites and goals of the study
Located on the southwestern tip of Mallorca Island, in the Balearic archipelago, sa Dragonera island has a surface of 276 ha ( Figure 2). Sa Dragonera has steep slopes that end in marine cliffs up to 350 m a.s.l. We selected the close islet of es Pantaleu as a control area since the islet has historically been free of rats. Climate is typical Aleppo pines (Pinus halepensis) and shrubs (e.g. Pistacia lentiscus, Phillyrea spp., Salvia rosmarinus, Ephedra fragilis and Cneorum tricoccon)). The flora of sa Dragonera is formed by at least 481 species, 23 of which are endemic to the Balearic Islands. Due to its tiny surface, es Pantaleu has a simpler plant community, with a few Aleppo pines. Apart from vegetation, sa Dragonera holds a community of invertebrates, mainly formed by detritivore, omnivorous, herbivorous and carnivorous species of the orders Hymenoptera, Araneae, Isopoda, Coleoptera and Mollusca (Gastropoda), an endemic omnivorous lizard species, insectivorous and frugivorous passerine birds, ornithophagous falcons and seabirds (gulls and shearwaters). Rats in Mediterranean islands are known to be omnivorous, consuming many arthropods, gastropods and a wide variety of plants, both F I G U R E 2 (a) Vegetation cover types in the study area and locations of the study plots (90 × 90 m) used for vegetation spatial pattern analysis. The map in the upper left corner shows the location of the study area in the Balearic Islands (western Mediterranean). (b) Elevation (meters a.s.l.), (c) slope (%) and (d) orientation maps of sa Dragonera island, all three variables obtained from a 2 × 2 m resolution DEM. Map lines delineate study areas and do not necessarily depict accepted national boundaries. fruits and seeds, leaves and vegetative parts (Cassaing et al., 2007).
They can prey on burrow nesting seabirds (mainly their chicks) but except for small petrels, the seabirds are still present in islands with a historical presence of rats, as are lizards endemic to the Balearic archipelago (Pérez-Mellado et al., 2008;Ruffino et al., 2009).
Here, we studied several features of vegetation at the two study sites, sa Dragonera and es Pantaleu. Given the availability and resolution of data at several scales (satellite images, orthophoto images and field data), we analysed: (a) the influence of climate variables on NDVI dynamics (1984-2021); (b) the presence of regime shifts on these dynamics; (c) the land cover dynamics (1956-2020); (d) the vegetation cover and spatial pattern dynamics, focusing on the comparison before-after rat eradication (2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018); and (e) the dynamics in plant diversity . We explain below the methods we followed for each of the analyses; unless stated, all analyses were performed for the two islands.

| Selection of satellite images for NDVI
To create the seasonal NDVI time series, we first obtained a land cover map (SIOSE, 2014; https://www.siose.es/). Then, the limits of the different land cover types were revised, and adjusted when necessary, by photo interpretation of a higher spatial resolution (15 cm) orthophoto (PNOA, 2018; https://pnoa.ign.es/). Accordingly, the six types of land cover found in the study area were either single or composite covers in which one single cover occupies at least 70% of the enclosure (open shrubland, mixture of shrubs and conifers, coniferous forest, riparian grove and cropland) or non-vegetable covers (cliffs) (Figure 2).
For the calculation of the NDVI time series , Landsat image processing was performed using Google Earth Engine (GEE), which is a cloud computing platform with large computational capabilities, for example, the calculation of time-series changes with the algorithm LandTrendr (for this last tool, see next section below) (Gorelick et al., 2017;Kennedy et al., 2018). The GEE is considered the most up-to-date method to obtain reliable NDVI time-series data (Kennedy et al., 2018;Xiao et al., 2020). The method also allowed us to apply the tool on the whole surface of the vegetation cover of interest (i.e. open shrublands, see below) in sa Dragonera and es Pantaleu. The image collection needs to be consistent through time, and not include noise (i.e. atmosphere, clouds and shadows), sensor differences or other anomalies. These data have been atmospherically corrected, and include a cloud shadow, water and snow mask produced using CFMASK (Kennedy et al., 2018;Zhu et al., 2015).
We selected the Landsat 4-5 TM, Landsat 7 ETM and Landsat 8 OLI Collection 2 level-2 atmospherically corrected surface reflectance images for our analysis (Holden & Woodcock, 2016). We generated the normalized difference vegetation index (NDVI) for the selected images as a proxy of vegetation productivity for sa Dragonera and the control of es Pantaleu. An increase in vegetation productivity with time (i.e. larger NDVI) may imply either a 'greener' vegetation (with the same coverage) or an increase in coverage (with the same greenness). We selected NDVI data from spring (March, April and May) and summer (July and August). In Mediterranean ecosystems, spring is the peak of plant productivity, whereas summer is the season when plant productivity is lowest.

| NDVI dynamics and influence of climate
We first built, for each season and site, generalized linear mixed models (GLMMs, Gaussian family), in which NDVI data were used as the dependent variable, several climate data were used as potential explanatory variables and the date of the image was introduced as a random effect. NDVI data were normally distributed for each season and site (Shapiro-Wilk tests not significant, results not shown).
A suitable model would allow us to assess the effect of climate on NDVI dynamics.
We tested the occurrence of a temporal trend in NDVI by adding the year as a continuous variable. Climate variables were obtained from a high-resolution Iberian dataset with a spatial resolution of 1 km and a weekly temporal resolution: evapotranspiration (ETO), relative humidity (RH), insolation (INS), mean temperature (Tmean), rainfall (RF) and that of that season the year before (RF_1). We also used the standardized precipitation-evapotranspiration index (SPEI), which is a multi-scalar index to identify drought duration and mag- March) (https://www.ncdc.noaa.gov), since NAO integrates several climate variables, which correlate with several processes in terrestrial (e.g. plant phenology) and marine ecosystems.
Since climate variables could show strong correlations, we made a selection based on PCA before using all of them in the GLMMs.
We checked the contributions (as the proportion of variance and eigenvectors) of each climate variable in each principal component (PC). We selected how many PCs were needed to describe a high cumulative percentage (>90%) of variance. In each of these PCs, we chose the variable with the highest score (irrespective of its positive or negative sign). Since PCs are orthogonal in the PCA, selected variables will be completely independent (non-correlated).
Model selection was made using the values of the corrected Akaike information criterion for small samples (AICc) (Burnham & Anderson, 2002). We also assessed the accuracy of the selected models using three methods: simple correlation analysis, fit of predicted values and mean absolute percentage error index (see details in section 'Accuracy of NDVI models' in Supplementary Material).
To calculate the temporal synchrony of NDVI in open shrublands between the two sites (sa Dragonera with rat eradication and the control site in es Pantaleu without rats), we calculated the proportion of local minima/maxima common to both time series and computed its statistical significance via Monte Carlo randomizations (Buonaccorsi et al., 2001;Gouhier et al., 2010).
To assess whether and when transitions in vegetation dynamics occurred, we applied the Landsat-based Detection of Trends in Disturbance and Recovery (LandTrendr) algorithm in GEE Kennedy et al., 2010), which constructs Landsat image time-series stacks (LTS) to detect temporal changes without the need to download and process the images (Kennedy et al., 2018).
LandTrendr generates a trend-fitting NDVI time series that remove year-to-year noise and test for transient dynamics in our study. We set the standard parameter (maximum number of segments = 6), which is considered suitable to detected short-and long-term changes (Kennedy et al., 2018). However, LandTrendr is designed to detect large changes resulting from abrupt environmental changes (e.g. fires, infrastructure building), but it may not identify more gradual changes especially when there is a temporal trend in the time series, as we found in our study (see Section 3). Thus, we also used a method designed to detect such type of ecological changes in a time series showing a trend. NDVI data were subjected to a 'dynamic shift detector' DSD tool, which uses an algorithm that identifies shifts in the time series with nonlinear dynamics (Bahlai & Zipkin, 2020). The algorithm implies an iterative approach to fitting subsets of timeseries data; later, it ranks the fit of breakpoint combinations using model selection based on the AIC values of the models, assigning a relative weight to each break. Finally, the method provides the rate of increase in each segment of the time series, in our case showing how NDVI would increase or decrease before and after each breaking point.

| Land cover dynamics
We elaborated a land cover mapping for the years 1956, 1984 and 2020 by photo interpretation, using the same six classes of land To analyse the dynamics of the land cover types, we carried out a cross-tab analysis to calculate an agreement index (Kappa index, KI, expressed as a percentage of change) for each matrix. A KI value = 0 indicates that the level of match is equal to that produced by chance, while a KI value = 100 indicates a complete match (i.e. no change) in land cover between two dates. Owing to the dynamics of perturbations at Dragonera (Figure 1), with a cessation of all harvesting in the mid-1970s, we expected a higher KI over the first period  than over the second (1984-2020).

| Vegetation cover and spatial pattern dynamics
We focused here on the differences between before and after rat eradication at sa Dragonera in 2011. For a balanced comparison, we obtained six high spatial resolution orthophotos of the study area The dynamics in plant cover were analysed by fitting a GLMM with a binomial distribution of errors. Specifically, year and topographic slope, as well as their interaction, were included as fixed terms in the model. The slope was introduced to take into account the spatial heterogeneity within sa Dragonera, since there is an important gradient of steepness (Figure 2c). The study plot was included as a random factor to deal with repeated measures. Similarly, the significance of differences in plant spatial structure PLR metrics over time was analysed by fitting separate linear mixed models (LMMs). Plant cover, year and topographic slope, as well as the interaction term between year and slope, were included as fixed factors in the models, while the study plot was included as a random factor. PLR data were arcsine-root transformed to meet normality assumptions.

| Temporal dynamics in plant diversity
We compiled floral taxon field inventory data collected with a similar methodology only at sa Dragonera in books and technical reports to assess how the number of taxa has changed over a century (1921-2021) (Oro et al., 2021;Sáez et al., 2013). Inventories were obtained in 13 different years, mostly (70%) since 1988. Since 2002, a standard methodology was set: the surface of the island was divided into 23 UTM 500 m 2 squares and each square was visited each sampled year at least once during several journeys (range 9-12) from March to October (Bibiloni et al., 2007). Cliffs were sampled from the distance using binoculars. Since data were sparse, especially until the late 1970s, we only performed a regression analysis to fit a logistic curve of the last 12 years with data, which would test a transient to a new stable state after cessation of anthropogenic perturbations. Most SPEI indexes (each with a different time lag) were correlated, especially those close in time. The SPEI that was most correlated with NDVI was the one with 12 months lag both in spring and summer (see Figure S4). Thus, only SPEI_12 was included in the PCA on climate variables. This analysis showed the first five PCs had an accumulated percentage of the variance of >95% and corresponded to the following climate variables: Tmean, NAOw, RF, RF_1 and RH for spring and Tmean, NAOw, RF, RF_1 and SPEI_12 for summer (see Figure S5 and Tables S1 and S2).

| NDVI dynamics and the influence of climate
There were only two land cover types present in the control islet of es Pantaleu: shrublands and rocky cliffs. Since the last cover is not suitable for studying NDVI dynamics, we focused our comparative analysis on the former cover. NDVI temporal dynamics of the differ-  Tables S3 and S4). Furthermore, the rainfall of the previous season (in spring for both islands) and the mean temperature (in summer for both islands) were the climate drivers that also play a statistically significant role in NDVI variability (see Tables S3 and S4). The temporal trends seem to stabilize in the last years, especially in sa Dragonera (Figure 4). The accuracy of the selected models for each island and season was very high for all accuracy indicators (see Table S5).
The two sites showed high synchronicity both in spring and in summer (proportion of common peaks κ = .7619, p = .010; κ = .8095, p = .010 respectively).
The increase in NDVI suggested nonlinear transient dynamics towards a new equilibrium (Figures 4 and 5). As we expected for a times series with an increasing trend, LandTrendr did not detect statistically significant segmentation neither using the standard parameter of the algorithm ( Figure 5) nor with other set of parameters (results not shown). On the contrary, DSD method detected two regime shifts occurring in spring in sa Dragonera, with an acceleration of NDVI that lasted a few years and a later deceleration that suggested stabilization of NDVI dynamics in the last years (Figures 4a and 5a) (see Tables S6 and S7). We detected only one F I G U R E 3 Mean temperature (Tmean, in °C), insolation (INS, in hours/day), evapotranspiration (ETO, in mm), accumulated rainfall (RF, in mm), relative humidity (RH, in %) and drought index SPEI for the season of the year before (SPEI_12) for spring and summer in the study site. A loess smoothing function (black line with standard error shaded in yellow) is fitted to each time series, and the value of r 2 for a linear model testing a temporal trend is also shown for each variable but SPEI_12 due to missing values. regime shift in es Pantaleu coinciding with the first one occurring in sa Dragonera, and there was an acceleration of NDVI in the last 20 years (Figures 4c and 5c). In summer, the NDVI showed a single regime shift in Dragonera, with also an acceleration in the last 15 years (Figures 4b and 5b). The NDVI of Pantaleu in summer showed similar dynamics than in Dragonera, with a synchronous regime shift in 2005. During this season, there was no indication of NDVI stabilizing in the last decades, but rather an acceleration in the increase in vegetation productivity (Figure 5b,d) (see Tables S8 and S9). Regime shifts were mostly coincident in the two islands and seasons and occurred earlier than rat eradication, which seems did not play a role in those shifts.

| Land cover dynamics
As expected, the coincidence (Kappa Index) in both vegetation surface by cover and location of the surfaces decreased between the first  and the second period (1984-2020) (from 73% to 69%), and it was lower when we considered the whole period  (KI = 52%). In other words, vegetation covers changed less over the first 29 years (most of this period was intensely harvested) than over the last 36 years when harvesting was already ended. An interesting result was that in 1956 when timber harvesting was intense (Figure 1), no coniferous forest cover was detected ( Figure S6). Over the first  and second  periods, part of the main vegetation covers (11% and 13% of the surface respectively) showed a F I G U R E 4 Temporal dynamics of NDVI  in open shrublands at sa Dragonera (panels a and b for spring and summer respectively) and the control islet of es Pantaleu (panels c and d for spring and summer respectively). We show a loess smoothing function to each time series to ease visualization of the trends. Arrows point to the last year with the presence of rats (and mice and rabbits) in sa Dragonera.  for open shrublands in sa Dragonera (red lines, a and b for spring and summer respectively) and the control site in es Pantaleu (green lines, c and d for spring and summer respectively). Arrows point to the last year with the presence of rats (and mice and rabbits) in sa Dragonera. Dashed lines show the years when a regime shift occurred; r shows the rate of increase in each segment of the NDVI time series. Orange lines show the trend-fitting NDVI time series using LandTrendr. similar shift to denser communities: open shrubland became a mixture of shrubs and conifers and this last cover partly became coniferous forests (Tables S10 and S11 and Figures S7 and S8).

| Vegetation cover and spatial pattern dynamics
First, we found a statistically significant positive relationship between NDVI and plant cover in open shrubland vegetation cover type in sa Dragonera (F 1,19 = 25.94; p < .001) ( Figure S9).
The interaction between year and slope was statistically significant for explaining changes in plant cover of the open shrubland cover in sa Dragonera ( 2 1 = 7.32; p < .01). We found that plant cover increased in study plots over time, but this increment was particularly large in lower slopes ( Figure S10). Contrarily, we found that the PLR index did not change significantly over time, and neither changed significantly with the slope nor with plant cover (p > 0.5). Similarly, the mark correlation functions for all years and study sites were lower than expected by the null model at r < 5 m (Figures S13 and S14). Again, these findings were similar for all study years and es Pantaleu islet (Figures S12 and S14).

| Dynamics of plant diversity
Despite the limited number of sampled years (13) during  and that records of plant species richness (as the number of plant taxa recorded) may be sensitive to annual both climatic conditions and sampling efforts, results show that there was a 3.5-fold increase in the number of species since the mid-70s ( Figure 6). That increase coincided with the end of harvesting practices, the removal of goats and the setting of some management actions over the last decades (e.g. sowing and cultivation of cereals in the old terraces). Besides, this number stabilized in the mid-90s and there was not a neat effect of rat eradication on the diversity of taxa recorded (Figure 6). Considering the whole time series but the single data from 1920, a logistic function fitted with this temporal dynamics (R 2 = .673, F 1,10 = 20.595, p < .001).
The increasing taxon diversity does not necessarily mean successive additions of taxa colonizing the island, but the net balance between those and the species that have likely disappeared (mainly annual plants from open habitats).

| The end of long-term harvesting by humans for vegetation dynamics
Small island ecosystems are suitable scenarios to study the influence of climate and perturbations on the dynamics of communities, which are composed of fewer species embedded in simpler, although still complex, food webs (Lopes et al., 2019;Markwell & Daugherty, 2002). We studied here the recent dynamical responses of several features of vegetation after the cessation of four millennia of intense anthropogenic impacts in sa Dragonera, a small island in the Mediterranean region. If we assumed that the historical information recorded in the close island of Mallorca applies to sa Dragonera, impacts started soon after the arrival of humans. Those impacts generated a loss of flora and fauna (e.g. the endemic both cave goat and giant dormouse) mainly by direct hunting, timber logging and the introduction of herbivores and rodents, and it likely triggered a transition to a new asymptotic equilibrium (Alcover et al., 2000;Wood et al., 2017). Historical chronicles show that the logging of pine forests by the inhabitants of the island was intense at least since the Middle Ages (CIM, 1996). Logging increased since the XVIII century when part of the timber was exported to Mallorca and the surface covered by forest suffered oscillations from overharvesting to recovery (Rosselló & Bover, 1980), typical of slow-fast dynamics in long transient phenomena (Hastings, 2004;Hastings et al., 2018). Botanists described the island as treeless in the 1850s and the 1900s; after that, there was a recovery of conifer F I G U R E 6 Dynamics of plant species richness in sa Dragonera during 1921-2021. Arrows point to the end of harvesting (timber logging and goat herbivory) and rat eradication. A loess regression line with standard errors is shown for the time series data since 1977.
forests until the 1940s, when most of them were cut down, and later recovery was boosted by pine reforestation in the 1970s when the reforestation area (ca. 40 ha) almost equalled the surface occupied by grown forests (ca. 45 ha) (CIM, 1996). This cyclic succession of vegetation dynamics caused by anthropogenic perturbation is similar to what has been observed in several ecosystems subjected to vegetation harvesting and management (Foley et al., 1996;Marrs & Hicks, 1986).
We found increasing trends in the dynamics of productivity, selforganization, fractal structure in the spatial patterns, forested cover and plant diversity, which were indicators of fast timescale dynamics occurring over the transient to an asymptotic stable state. For instance, there was a predominant spatial pattern of aggregation in small patches showing facilitation rather than competition in plant communities, which would occur after a perturbation, far from density dependence and the carrying capacity of the system. Sparse data from 1921 and 1956 confirmed that, before the end of intense harvesting, vegetation showed an ecological state of a perturbed system, lacking forested patches and with 3.5-fold of lower plant diversity (see also Picornell-Gelabert et al., 2021). Vegetation cover increased more in flat than in step areas likely because the former were more impacted by anthropogenic activities and hold higher rat densities (Mayol et al., 2012). Spatial heterogeneity in the effects of rats has been found in a range of ecosystems, including other islands (Himsworth et al., 2013;Igual et al., 2006;Rayner et al., 2007).

Strikingly, high synchronicity in vegetation productivity between sa
Dragonera and the control site (free of rats) showed that rat eradication did not alter those transient dynamics. A previous study, which analysed NDVI dynamics in sa Dragonera and used a control site with rats and goats (Ibañez-Álvarez et al., 2022), showed also synchronicity between the sites. Altogether, these results suggest that neither rats nor goats were the main drivers of the transient vegetation dynamics we observed, but rather the synergetic effects of timber logging and natural-and human-induced fires (the last intended to favour domestic herbivore grazing) that characterized the harvesting practices before the abandonment of a rural-dominated economy (Burjachs et al., 2017;Mus, 1995). In the Mediterranean region, transient vegetation dynamics have been more related to anthropogenic impacts and societal collapses than to climate variability, especially in the last four millennia (Carrión et al., 2001;Fyfe et al., 2019;Pausas, 1999).

| An evolutionary view of vegetation dynamics
Vegetation dynamics in sa Dragonera show similar increasing trends to those recorded in many regions of the industrialized world due mainly to climate warming, land-use changes and rural abandonment (Cervera et al., 2019;García et al., 2019;Scheffers et al., 2016;Winkler et al., 2021). In islands where conservation actions have eradicated invasive herbivore mammals, the increase in plant communities and forested cover have been even larger (Kappes et al., 2021;Luna-Mendoza et al., 2019;Médail, 2017). A singularity of sa Dragonera was that owing to its small size, human settlements were historically nil or much reduced, which likely caused the lack of introduced carnivores (Bover & Alcover, 2008) and the persistence of endemic lizards, which play an important role as shrub seed dispersers (Traveset & Riera, 2005). Over more than five million years, and before the arrival of humans, ecological communities and their interactions evolved under climate changes and environmental stochasticity throughout the Pliocene, the Pleistocene and the Holocene (e.g. glaciations, inter-glacial periods and the Holocene warm Optimum) (Alcover et al., 2000;Burjachs et al., 2017)

| Ecosystem dynamics after the end of harvesting
Having in mind the evolutionary singularities of sa Dragonera, we may wonder whether the transient vegetation dynamics of sa Dragonera are now on a path returning to the asymptotic equilibrium before the arrival of humans. Since there is no record of fossil fauna in sa Dragonera during the Pliocene, we cannot rule out that populations of rodents and bovid suffered extinction-colonization processes, as it was recorded for larger islands in the Balearic archipelago (Mas et al., 2018). Despite this lack of detailed information, we hypothesize that the ecosystem had a transition from predatorregulated populations of an endemic herbivore and an omnivore to populations of domestic goats and rats likely attaining higher population densities than the native ecological counterparts by the lack of predators and, for rats, the exploitation of anthropogenic food subsidies (Oro et al., 2013;Ruffino et al., 2013). In the last four millennia, environmental stochasticity, competition, harvesting rates and density-dependent processes related to vegetation dynamics shaped the population dynamics of rats and goats. By being omnivores, rats had not only direct but also indirect effects on vegetation dynamics: they likely regulated the density and cycles of phytophagous insect outbreaks (e.g. bark beetles) and indirectly the dynamics of seed dispersal by predation of frugivorous passerines, which are now released from these predator-prey dynamics. Qualitative data on plant composition recorded at sa Dragonera in the last decade show that the density of saplings of several shrubs has increased (e.g. mallow Lavatera spp., dwarf palm), whereas annual plants and grasslands of open habitats are decreasing and even disappearing (Oro et al., 2021. New environmental policies have also affected other components of the ecosystem shaping vegetation dynamics. For instance, ground-breeding gulls attained large population densities by being released from historical harvesting of their eggs and by anthropogenic food subsidies (e.g. dumps, fishing discards) in the last four decades (Payo-Payo et al., 2015). As a consequence, there was increased vegetation growth by a higher supply of nutrients and larger seed dispersal (i.e. colonization processes from the outside) Calviño-Cancela, 2011;Jones, 2010;Rosselló & Bover, 1980;Ruffino et al., 2013). The nutrient supply by gulls drastically dropped since 2010, when the dump where gulls mostly fed was closed and halved their population density (Payo-Payo et al., 2015). That supply was counterbalanced by the recolonization of the island by burrow-breeding seabirds decimated by rats, such as shearwaters (Martin et al., 2000;Morgan et al., 2013). Despite that assessing the net balance between all interactions in the food web of sa Dragonera for vegetation dynamics is challenging, it seems clear that larger forested areas, denser shrublands and loss of open habitats would characterize the new ecosystem state.

| Ecological stability of vegetation after island ecosystem restoration
Socio-economic dynamics of land use and conservation management eradicating alien species (mainly domestic ungulates, carnivores and rodents) have released many small island ecosystems from most anthropogenic impacts. This represents an opportunity to study ecosystem stability and transient dynamics after long-term perturbations. After five decades of lack of historical harvesting and free of rats, goats, mice and rabbits, vegetation communities at sa Dragonera would be in the final stage of a transient towards a stable long-term equilibrium. The eventual asymptotic dynamics of vegetation would be regulated by fewer ecological interactions (lacking large herbivores and omnivores), climate variability (mainly the severity and occurrence of droughts), natural fires and density-dependent processes set by the carrying capacity of the island (i.e. availability of soil, nutrients and water). When we frame the vegetation dynamics in terms of stability and dynamical systems (Hastings, 2004, p. 200;Holling, 1973, p. 197), the regime shifts we found show that the vegetation state is moving from a domain of attraction characterized by open shrubland to a state dominated by denser shrubs and pine forests (see a similar dynamics in Glendening, 1952). The stability of the new equilibrium will be enhanced by two factors: first, the adaptation of Mediterranean vegetation communities (from annual plants to shrubs and trees) to be resilient to fires, droughts and strong climate seasonality, including extreme climate events (Carrión et al., 2010;Pausas, 1999). Second, and perhaps more importantly, conservation actions target the avoidance of both rat re-invasion and occurrence of fires and the control of insect outbreaks affecting shrubs and trees. In that scenario, anticipating the dynamics of populations of the most valuable taxa for conservation, such as endemic lizards, birds, invertebrates and plants is challenging. Previous studies found unforeseen effects of management actions, including the eradication of alien species (Ausden et al., 2001;Martínez-Abraín et al., 2004). Return to the original domain of attraction set by open shrubland and the co-evolution of communities including herbivores and omnivores, can only be made by an explicit reduction of the trees and shrubs. This goal, which should increase ecosystem resilience, may be attained by grazing goats at low densities, vegetation clear-cut management and controlled fires (Mancilla-Leytón et al., 2013;Marull et al., 2015;Naveh, 1990).

CO N FLI C T O F I NTE R E S T
The authors declare no competing interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are openly available in CSIC Digital Repository (DIGITAL.CSIC) at http://hdl.handle. net/10261/ 266484.