Satellite-Based Monitoring of Primary Production in a Mediterranean Islet Post Black Rat Eradication

: Invasive rodents have a detrimental impact on terrestrial ecosystem functioning, this is often exacerbated on small islands. Rat eradication campaigns are often used to deal with this environmental perturbation given their classiﬁcation as invasive species. Studies assessing the effects of rodent control at ecosystem scale are scarce and thus little is known about the subsequent functional response of vegetation subsequent to rat control. In this work, we use remote sensing to assess the effects of black rat ( Rattus rattus ) eradication on Mediterranean vegetation productivity in the Sa Dragonera Islet, Mallorca (Spain). Rats feed on seeds, sprouts, and leaves of woody vegetation and hence we expect primary production to increase nine years after the rodenticide campaign. The Break Detection approach for additive season and trend (BFAST method) was adopted to examine changes in vegetation density before and after the eradication campaign in Sa Dragonera Islet (Balearic Islands), using a temporal series of monthly NDVI data extracted from Landsat imagery. The same temporal trends were examined for a control zone where no rat eradication took place, in order to control for weather-driven changes. The results of this study revealed changes across the 21-year monthly NDVI time series. However, the dates, magnitude, and trend of these changes could not be explicitly attributed to the action of rats, when compared to the historical changes on the islet and the changes found to co-occur within the control zone. These ﬁnding could, perhaps, be explained by the high resilience of Mediterranean shrubs to browsing including that of rat invasion. However, the results from the study appear to show that rat damage on speciﬁc plant species, with little contribution to global NDVI values, would be overshadowed by the effects of broader environmental factors in this remote sensing approach. The results suggest that the current passive restoration scheme imposed following eradication is not sufﬁcient for effective ecosystem restoration. The BFAST results reveal that the primary production in the treatment area is sensitive to water cycles but not immediately to cessation of rat activity. In the short term there was no abrupt increase in vegetation primary productivity as a response to the cessation of rat browsing with this signiﬁcant ecological event being overshadowed by the state of water emergency that occurred close to the time of the eradication campaign. Nor were changes observed in the long term, since facing the aforementioned environmental stress events, the recovery rates in the treatment zone did not overcome those found in the control zone, rather, they continued to show similar trends. Changes may have only occurred at a small spatial scale, as evidenced by the reported ﬁeld observations of the anecdotal studies in the area. Thus, given the spatial resolution of the core satellite data these changes had little impact on the averaged NDVI produced over the islet. Unlike other studies, changes were not widespread over the island’s surface, so it is not possible to detect them at the scale at which Landsat sensors record. We suggest that the non-response of NDVI to rat eradication on Dragonera Islet is due to a high resistance of vegetation to rat predation achieved during a long period of colonization in an already naturally stressed environment.


Introduction
Invasive alien species (IAS) are those that have reached new geographic areas by way of human introduction, and are populations which can survive, reproduce and spread over the natural environment leading to major impacts on the environment or society [1]. IAS can cause serious disturbance to other indigenous species and the ecosystems they invade, often altering ecosystem structure and function, trophic relationships (e.g., plant-animal interactions), and reducing the biodiversity of the invaded ecosystems [2][3][4][5].
Island communities which have evolved in isolation are typically the most vulnerable to IAS [6]. One of the most harmful invasive species are rodents belonging to the genus Rattus, in particular, the black rat (Rattus rattus), brown rat (R. norvegicus) and Pacific rat (R. exulans) [7][8][9]. The impact of rats is often density-dependent [10][11][12] and regulated by primary production [13][14][15], thus, in many biomes rat populations fluctuate seasonally [16]. Records of black rat density are highly variable ranging from 36.4 rats ha −1 in New Zealand islands [15], to 50 rats ha −1 in some Mediterranean islands [17]. In the Mediterranean basin the black rat is the most invasive species due to its broad tolerance to dry climates [18,19]. These species have already invaded 80% of the world's islands [12,[20][21][22] and they cause a cascade of damage to ecosystems, resulting in the decline and extinction of native birds, mammals, reptiles, invertebrates and plants of many archipelago ecosystems [7,8,18]. The black rat tends to consume more plants than animals [23], showing a preference for seeds, sprouts, and leaves of adult plants [24,25], and thus directly affecting photosynthetic and reproductive parts. As a result, black rats affect the composition and structure of plant communities, plant recruitment, and the viability of the plant population on the islets and islands they infest. This has been recorded in some islands of New Zealand [26,27], and in islands and small islets of Spanish archipelagos (Canarias and Balearic Islands) where the presence of black rat has altered the structure of the vegetation; depressing or limiting some plant species and favoring others [19].
In islands where rats are removed at an early stage of their invasion, vegetation recovers quickly [28,29]. However, in places with prolonged history of rat colonization, recovery is much slower due to critical changes in vegetation structure [30,31]. Rodent eradication is a common measure [8] for recovering invaded areas [32] with accompanying environmental monitoring necessary to evaluate the effectiveness of rat eradication campaigns in terms of ecosystem restoration [33]. However, ecosystem recovery is often affected by other confounding factors (e.g., climatic events) [34,35], or is limited to specific plant or animal species [36,37]. Conversely, the lack of environmental information prior to rat invasion [38] can hamper the evaluation of ecosystem restoration programs [39]. In fact, ecosystem monitoring post rodent eradication is uncommon [13,40]. Remote sensing therefore has a clear role to play in assessing ecosystem restoration post rodent eradication [41] given its spatial and temporal extents.
Satellite sensors are capable of providing land cover information at different spatial, spectral, and temporal resolutions, allowing the monitoring of vegetation extent and health, and through this indicate areas of recovery post natural or human-induced disturbance [42,43]. A range of spectral vegetation indices, such as the Normalized Difference Vegetation Index (NDVI), can be used for these purposes [44]. NDVI is a proxy for photosynthetic activity [45], vegetation productivity, aboveground biomass [46], and vegetation dynamics [47,48], and has been shown to reliably capture ecological response to environmental change including human and animal-driven land degradation [44]. In most terrestrial biomes, long-term NDVI time series follow a non-linear pattern that comprises alternating cycles of greening (increasing NDVI) or browning (decreasing NDVI). The dating and quantification of frequency and magnitude of these NDVI cycles is generally used as a proxy for ecosystem dynamics [44,48,49]. NDVI time series are particularly useful in periods longer than 10 years [50], as this makes it feasible to analyze the seasonal and trend components over time. The nature of the change in NDVI time series can be interpreted according to the affected component. For example, modifications of the seasonal oscillation could indicate phenological changes, while changes in the trend, including the magnitude, may point to a disturbance event [43], precipitation variability [51][52][53] or drought event [53,54].
In this work, we evaluate the impact of a rodent eradication campaign in a Mediterranean scrubland located in the Baleares Archipelago (western Mediterranean basin, Spain) using twenty-one years of monthly NDVI data for the Sa Dragonera and for a nearby control zone on the island of Mallorca. Changes in vegetation productivity, before and after the eradication campaign, across the two zones was examined using the breaks for additive seasonal and trend analysis (BFAST, [43,55]). BFAST is especially useful here to address the lack of information regarding the magnitude of change on primary production after deratization management, because unlike other methods, BFAST does not restrict the data to a specific season since it considers a seasonal fit over the entire time series. Additionally, its change detection approach based on signal-to-noise ratio does not require setting a change threshold. BFAST has the sensitivity to differentiate the normal phenological cycle from abnormal changes and has been used to monitor the recovery of vegetation [56,57] after a broad range of environmental catastrophes including wildfires, insect outbreaks [58] and floods [59]. Due to the devastating effect that black rat infestations can have on small islet ecosystems reported worldwide [19] and the ability of BFAST to detect changes at the local scale [55,60,61], we hypothesized that primary productivity on Sa Dragonera Islet would differ before and after the rodenticide campaign; due to vegetation regeneration following black rat eradication. In particular, we expect to detect an increase (abrupt and/or gradual) in the greenness of the main vegetation due to the natural vegetation recovery following rodent eradication after controlling for climatic conditions.

Study Area
The focal geographic area of the study is in the Baleares Archipelago (western Mediterranean basin, Spain). It comprises Sa Dragonera Islet, referred to hereafter as treatment zone. In addition, a control zone was included, some 800m away on Mallorca Island, due to its similarity in orography, lithology and vegetation type. Importantly, both zones are an extension of the UNESCO designated World Heritage Site of the Serra de Tramuntana Mountain Range. The treatment zone is a 276 ha islet extending 350 m above sea level, with the north-west exposure terminating in a cliff face.
The climate at both sites is Mediterranean, with annual average temperatures above 16 degrees and little annual rainfall (<400 mm); but with high inter-annual variability. According to the Spanish Meteorological Agency (AEMET) [62], drought periods typically last 5 months; yet these events are partially mitigated by maritime influence. Limestone is the predominant bedrock and overlying soils are poorly developed or skeletal. The main vegetation consists of a mosaic of scrub communities. There are tree patches, and the herbaceous taxa are scarce (other than non-palatable plants such as Urginea maritima and Crithmum maritimum) due to the poor climatic and edaphic conditions. Sclerophyllous woody species of Mediterranean optimum predominate, such as Pinus halepensis, Pistacia lentiscus, Olea europaea, Anthyllis cytisoides, Chamaerops humilis, Cistus monspeliensis, Cneorum tricoccon, Ephedra fragilis, Euphorbia dendroides, Phillyrea angustifolia [63,64]. Most of these species can resprout following disturbance, exhibiting some resilience to fires or/and browsing, which are typical characteristics of sclerophyllous flora of the Mediterranean basin [65][66][67].
The islet is currently uninhabited; however, human activity has shaped its landscape as well as that of Mallorca Island over the last five millennia, fostering deforestation and soil erosion processes [63,68]. Natural resources, such as woodcutting, crop growing, and domestic animal breeding, were common until the 1975 [69,70]. Feral goats (Capra hircus) were eradicated in the 1980s, European wild rabbits (Oryctolagus cuniculus), and house mouse (Mus musculus) were also present in the treatment zone, yet black rats have been by far the most common invasive mammal in the islet reaching up to 50 individuals ha −1 [71]. This density of rats equals the maximum reported in Mediterranean islands [17,72] and is higher than that historically recorded in the islands of New Zealand (36.4 rats ha −1 ) [15].
Sa Dragonera is the main islet of an archipelago which was awarded the status of Natural Park in 1995 and forms part of the Natura 2000 Special Protection Area and Site of Community Importance (SPA and SCI, EU Birds and Habitats Directive) [73]. The majority of the islet is considered a reserve area devoted to the conservation of natural values in the absence of active management measures [74]. Field observations revealed severe damage by rodents on vegetation cover, for example, in 2010 partial or total damage on branches of juvenile individuals of Ficus carica, Ceratonia siliqua and Olea europaea var. sylvestris were reported ( Figure 1). Likewise, the endangered populations of Balearic shearwater (Puffinus mauretanicus) were impacted by rat predation [71]. Hence, conservationists have performed several ground-based rodent eradication attempts between 2001 and 2008 using rodenticide baits in the most accessible areas, but with little long-term success as rat populations have recovered to high densities. In February 2011 the local government organized an eradication campaign using Brodifacoum rodenticide. Brodifacoum baits were aerially distributed twice over the entire islet at rates of 14 kg ha −1 [71]. Since this eradication effort, the treatment zone has been rodent free [75,76], despite the high risk of reinvasion due to the proximity of the coast of Mallorca [77]. Four years after Brodifacoum baiting, unsystematic and occasional field observations have reported the increase of vegetation cover [78], however evidence of this recovery remains lacking.
Natural Park in 1995 and forms part of the Natura 2000 Special Protection Area and Site of Community Importance (SPA and SCI, EU Birds and Habitats Directive) [73]. The majority of the islet is considered a reserve area devoted to the conservation of natural values in the absence of active management measures [74]. Field observations revealed severe damage by rodents on vegetation cover, for example, in 2010 partial or total damage on branches of juvenile individuals of Ficus carica, Ceratonia siliqua and Olea europaea var. sylvestris were reported ( Figure 1). Likewise, the endangered populations of Balearic shearwater (Puffinus mauretanicus) were impacted by rat predation [71]. Hence, conservationists have performed several ground-based rodent eradication attempts between 2001 and 2008 using rodenticide baits in the most accessible areas, but with little long-term success as rat populations have recovered to high densities. In February 2011 the local government organized an eradication campaign using Brodifacoum rodenticide. Brodifacoum baits were aerially distributed twice over the entire islet at rates of 14 kg ha −1 [71]. Since this eradication effort, the treatment zone has been rodent free [75,76], despite the high risk of reinvasion due to the proximity of the coast of Mallorca [77]. Four years after Brodifacoum baiting, unsystematic and occasional field observations have reported the increase of vegetation cover [78], however evidence of this recovery remains lacking.
The Sierra Tramuntana in Mallorca is also a landscape shaped by human activity, which suddenly ceased from the 1960s in the rural mountain areas, becoming focused on rapid development of the tourism sector in the urbanized areas. The abandonment of the rural sector led to the cessation of domestic animal grazing but also to the proliferation of the feral goat which continues today [79,80]. The Sierra Tramuntana in Mallorca is also a landscape shaped by human activity, which suddenly ceased from the 1960s in the rural mountain areas, becoming focused on rapid development of the tourism sector in the urbanized areas. The abandonment of the rural sector led to the cessation of domestic animal grazing but also to the proliferation of the feral goat which continues today [79,80].

NDVI Time Series Construction
To identify whether the rodenticide eradication campaign affected primary production in the treatment zone, we compare its NDVI time series with that in the control zone. NDVI time series are here based on Landsat surface reflectance (SR) data for a period ranging from 1999 to 2020. February 2011 was set as the base date for detecting changes due to rat eradication. Therefore, the study period was divided into: (i) a historical period of 12 years where vegetation was impacted by rats (1999-2010) and (ii) a recovery period of 9 years post rat eradication (2012-2020).

Satellite Images
This study used the Landsat satellite archive due to the temporal continuity across the analysis time period. Landsat imagery Collection 1 provided by USGS (United States Geological Survey) was accessed via the web site https://earthexplorer.usgs.gov/ on 14 November 2020. Landsat Collection 1 images are calibrated and aligned across all Landsat satellites. 333 Landsat SR images were obtained for the period January 1999 to December 2019 from Landsat 5 TM, Landsat 7 ETM+ and Landsat 8 OLI sensors. SR images are a Landsat level-2 product that provide the bottom of atmosphere surface spectral reflectance as it would be measured at ground level (i.e., with atmospheric correction applied). The availability of a high number of images favors the temporal consistency of the radiometric data across time series, since gaps produced due to the presence of clouds or due to SLC off-Data error, inherent to Landsat 7 [81] acquisitions post May 2003, can be replaced with optimal values recorded by another sensor at a close date ( Figure 2).

NDVI Time Series Construction
To identify whether the rodenticide eradication campaign affected primary production in the treatment zone, we compare its NDVI time series with that in the control zone. NDVI time series are here based on Landsat surface reflectance (SR) data for a period ranging from 1999 to 2020. February 2011 was set as the base date for detecting changes due to rat eradication. Therefore, the study period was divided into: (i) a historical period of 12 years where vegetation was impacted by rats (1999-2010) and (ii) a recovery period of 9 years post rat eradication (2012-2020).

Satellite Images
This study used the Landsat satellite archive due to the temporal continuity across the analysis time period. Landsat imagery Collection 1 provided by USGS (United States Geological Survey) was accessed via the web site https://earthexplorer.usgs.gov/ on 14 November, 2020. Landsat Collection 1 images are calibrated and aligned across all Landsat satellites. 333 Landsat SR images were obtained for the period January 1999 to December 2019 from Landsat 5 TM, Landsat 7 ETM+ and Landsat 8 OLI sensors. SR images are a Landsat level-2 product that provide the bottom of atmosphere surface spectral reflectance as it would be measured at ground level (i.e., with atmospheric correction applied). The availability of a high number of images favors the temporal consistency of the radiometric data across time series, since gaps produced due to the presence of clouds or due to SLC off-Data error, inherent to Landsat 7 [81] acquisitions post May 2003, can be replaced with optimal values recorded by another sensor at a close date ( Figure 2).
Once Landsat SR images were downloaded, R script was coded to search for the red and NIR spectral bands of each satellite observation, combining them according to the NDVI formula [82] using the "Raster" package [83] in the R statistical software [84]. Thus, NDVI layers were produced for all observations. All layers were stacked creating a multivariable raster dataset called RasterStack object.  Once Landsat SR images were downloaded, R script was coded to search for the red and NIR spectral bands of each satellite observation, combining them according to the NDVI formula [82] using the "Raster" package [83] in the R statistical software [84]. Thus, NDVI layers were produced for all observations. All layers were stacked creating a multi-variable raster dataset called RasterStack object.

Plot Selection Procedure
Time series are created by averaging NDVI values from selected sampling plots. The plots were located using a pair of UTM coordinates indicating the center and extended to the spatial resolution of a Landsat image (i.e., 900 m 2 ). To select sampling plots where vegetation shows similar temporal dynamics of NDVI, a sampling procedure was used to choose plots meeting the following features: (I) same vegetation cover, (scrubland, because it is the principal vegetation cover of the islet), based on Sentinel-2 Global Land Cover (S2GLC, [85]); (II) similar NDVI values avoiding scattered scrub areas and cliffs. For this purpose, a NDVI composition map was created where each pixel represented the 5-year average NDVI prior to the rodenticide event (i.e., from 2006 to 2011). This map was classified into the major categories (i.e., coniferous forest, scrubland, scattered thicket and bare soil) using the Jenks natural breaks classification method to maximise the variation between classes and to best group similar values [86]. With the help of a high resolution orthophoto (15 cm spatial resolution aerial orthophotography from the PNOA (National Plan for Aerial Orthophotography) available in http://www.scne.es/ #ORTOPNOA) on 17 October 2020, we photo-interpreted and established that the values of the target vegetation cover in the NDVI composition map ranges from 0.45 to 0.62. To address NDVI differences due to slope, altitude and orientation issues [87], sampling plots were selected if they were between 50 to 300 m asl, had a slope less than 30 degrees, and had an east or south facing-orientation. This selection was carried out using a 5 × 5 m 2 resolution Digital Elevation Model (DEM) raster layer in ArcGIS 10.7. A workflow of the sampling procedure is shown in Figure 3. A total of fifty-four sampling plots, 29 in the treatment zone and 25 in the control zone, met these criteria ( Figure 4).

Monthly Landsat NDVI Time Series Construction
For each sampling plot, NDVI values were extracted from a multilayer raster object (Section 2.2.1) using the "raster" 3.4-10 package [83] in the R statistical software [84]. The cleaning, gap-filling and smoothing procedure summarized below was followed: a) NDVI values outside the range of the target vegetation index range, below 0.2 or over 1, were removed. The NDVI values were then selected according to maximum-value composite (MVC) criterion. b) a gap-filling with linear interpolation between neighboring values [88] was performed, to fill NDVI gaps (31 observations in the treatment zone and 32 in Mallorca Island, the control zone). c) the Savitzky-Golay [89] smoothing algorithm was applied to smooth out the noise in both time series. It is appropriate to apply it on seasonal data [90] since it follows withinseason variations and therefore captures subtle dynamics during seasons [91].
By applying the steps described here, clean and consistent time series were constructed for each treatment to feed the BFAST algorithms ( Figure 5), since a temporal consistency of observations is critical for understanding ecosystems dynamics [92].

Monthly Landsat NDVI Time Series Construction
For each sampling plot, NDVI values were extracted from a multilayer raster object (Section 2.2.1) using the "raster" 3.4-10 package [83] in the R statistical software [84]. The cleaning, gap-filling and smoothing procedure summarized below was followed: (a) NDVI values outside the range of the target vegetation index range, below 0.2 or over 1, were removed. The NDVI values were then selected according to maximum-value composite (MVC) criterion.
(b) a gap-filling with linear interpolation between neighboring values [88] was performed, to fill NDVI gaps (31 observations in the treatment zone and 32 in Mallorca Island, the control zone).
(c) the Savitzky-Golay [89] smoothing algorithm was applied to smooth out the noise in both time series. It is appropriate to apply it on seasonal data [90] since it follows within-season variations and therefore captures subtle dynamics during seasons [91].
By applying the steps described here, clean and consistent time series were constructed for each treatment to feed the BFAST algorithms ( Figure 5), since a temporal consistency of observations is critical for understanding ecosystems dynamics [92].

Historical Time Series Analysis
To detect and characterize changes in the ecosystem dynamics, the BFAST protocol described in [43] was applied to the monthly NDVI time series. BFAST interactively decomposes time series into trend, seasonal, and remainder components [43,55]. The iteration begins using a Seasonal-Trend decomposition (STL) method. STL first detrends the time series by subtracting the trend from the observed values resulting in a new time series that exposes seasonality. The detrended time series is then used to compute the average seasonality. The amount of NDVI beyond the combination of the seasonal and trend components is considered random noise or the remainder component ( Figure 6). Subsequently, BFAST function detects breakpoints by applying ordinary least squares residuals-based moving sum (OLS-MOSUM). This tests the null hypothesis that the linear trend of the NDVI values remain constant over the observations against the alternative of significant changes (p < 0.05). The optimal number of breakpoints is determined by minimizing the Bayesian Information Criterion (BIC) and the position of the breakpoints (dates and its confidence intervals) are chosen by globally minimizing the residual sum of squares. BFAST iterates along these steps until the number and position of the breakpoints remains constant. The basis and the process carried out by BFAST is described in detail in [43] and is available in the BFAST package for R (https://cran.r-project.org/web/packages/bfast/index.html) on 9 September, 2021. All BFAST analyses were performed using R statistical software [84].

Historical Time Series Analysis
To detect and characterize changes in the ecosystem dynamics, the BFAST protocol described in [43] was applied to the monthly NDVI time series. BFAST interactively decomposes time series into trend, seasonal, and remainder components [43,55]. The iteration begins using a Seasonal-Trend decomposition (STL) method. STL first detrends the time series by subtracting the trend from the observed values resulting in a new time series that exposes seasonality. The detrended time series is then used to compute the average seasonality. The amount of NDVI beyond the combination of the seasonal and trend components is considered random noise or the remainder component ( Figure 6). Subsequently, BFAST function detects breakpoints by applying ordinary least squares residuals-based moving sum (OLS-MOSUM). This tests the null hypothesis that the linear trend of the NDVI values remain constant over the observations against the alternative of significant changes (p < 0.05). The optimal number of breakpoints is determined by minimizing the Bayesian Information Criterion (BIC) and the position of the breakpoints (dates and its confidence intervals) are chosen by globally minimizing the residual sum of squares. BFAST iterates along these steps until the number and position of the breakpoints remains constant. The basis and the process carried out by BFAST is described in detail in [43] and is available in the BFAST package for R (https://cran.r-project.org/web/packages/bfast/index.html) on 9 September 2021. All BFAST analyses were performed using R statistical software [84].
BFAST transforms the NDVI time series into a piecewise model, whose maximum number of linear regressions depends on the parameter h. The h parameter sets the minimal segment size between potentially detected breaks in the trend model given as a fraction relative to the sample size. The sample size of each of our time series consists of 252 observations (i.e., 21 years at a monthly frequency). We established the value of h = 1/7 in order to set significant trend periods following the recommendations of [59]. For our dataset, this implies that a maximum of seven linear regressions with 252/7 = 36 minimum observations each (i.e., a minimum of three years). The difference between the NDVI value at the end of a preceding linear model and the intercept of the following linear model defines the magnitude of the abrupt change and the slope of each linear model defines the gradual change between breakpoints. BFAST transforms the NDVI time series into a piecewise model, whose maximum number of linear regressions depends on the parameter h. The h parameter sets the minimal segment size between potentially detected breaks in the trend model given as a fraction relative to the sample size. The sample size of each of our time series consists of 252 observations (i.e., 21 years at a monthly frequency). We established the value of h = 1/7 in order to set significant trend periods following the recommendations of [59]. For our dataset, this implies that a maximum of seven linear regressions with 252/7 = 36 minimum observations each (i.e., a minimum of three years). The difference between the NDVI value at the end of a preceding linear model and the intercept of the following linear model defines the magnitude of the abrupt change and the slope of each linear model defines the gradual change between breakpoints.
For better understanding of the vegetation dynamics and to avoid misinterpretation of the effect of rat removal, we compared the BFAST outcomes between zones and with historical environmental records. Because changes in hydrological regimes are considered the main drivers of vegetation dynamic changes [52,54], and closely related to shifts in NDVI time series in water limited ecosystems [93,94], we compared the BFAST results with the hydrological drought index (HDI). This is reported by the General Directorate of Water Resources (DGRH) of the government of the Balearic Islands [95] in the Hydrological Demand Unit (HDU) of Sierra Tramuntana Sur. This corresponds to our study area. In the report they categorize HDI in 4 groups which describe the state of water resources. These categories, from highest to lowest hydric abundance are: normality, pre-alert, alert and emergency. The years 1999, 2001, 2004, 2007 and 2011 are key dates to consider due to the area entering alert or emergency states. Similar to the majority of zones in Mallorca, as a result of the high rainfall in 2009 and 2010, our study area was at optimal HDI levels, although from 2015 onwards a rapid decline was observed, leading to the alert situation in late 2015.
The BFAST outputs are compared between zones. This, on the one hand, gives us an idea of whether the process of selection of plots was efficient at the time of selecting areas whose vegetation follows the same dynamics. On the other hand, the comparison of the dates when the breaking points occur tries to find the singularity of the effect of the eradication campaign. In addition, the magnitude of the breakpoints is used as a proxy for the vegetation resistance facing a disturbance event; and the slope of the linear regressions as For better understanding of the vegetation dynamics and to avoid misinterpretation of the effect of rat removal, we compared the BFAST outcomes between zones and with historical environmental records. Because changes in hydrological regimes are considered the main drivers of vegetation dynamic changes [52,54], and closely related to shifts in NDVI time series in water limited ecosystems [93,94], we compared the BFAST results with the hydrological drought index (HDI). This is reported by the General Directorate of Water Resources (DGRH) of the government of the Balearic Islands [95] in the Hydrological Demand Unit (HDU) of Sierra Tramuntana Sur. This corresponds to our study area. In the report they categorize HDI in 4 groups which describe the state of water resources. These categories, from highest to lowest hydric abundance are: normality, pre-alert, alert and emergency. The years 1999The years , 2001The years , 2004The years , 2007The years and 2011 are key dates to consider due to the area entering alert or emergency states. Similar to the majority of zones in Mallorca, as a result of the high rainfall in 2009 and 2010, our study area was at optimal HDI levels, although from 2015 onwards a rapid decline was observed, leading to the alert situation in late 2015.
The BFAST outputs are compared between zones. This, on the one hand, gives us an idea of whether the process of selection of plots was efficient at the time of selecting areas whose vegetation follows the same dynamics. On the other hand, the comparison of the dates when the breaking points occur tries to find the singularity of the effect of the eradication campaign. In addition, the magnitude of the breakpoints is used as a proxy for the vegetation resistance facing a disturbance event; and the slope of the linear regressions as a proxy of the recovery rate (greening when positive rate or browning when negative rate) [96].

Historical Time Series Results
BFAST decomposed the original NDVI data (Yt) into seasonal (St), trend (Tt) and remainder (et) components (Figure 7). The algorithm identified four breakpoints at both the treatment and control zones in the NDVI trend component (Figure 7a,b, respectively), two of them following the eradication campaign. The seasonal component showed no abrupt changes, and the seasonal amplitude was approximately 0.05 NDVI.

Historical Time Series Results
BFAST decomposed the original NDVI data (Yt) into seasonal (St), trend (Tt) and remainder (et) components (Figure 7). The algorithm identified four breakpoints at both the treatment and control zones in the NDVI trend component (Figure 7a,b, respectively), two of them following the eradication campaign. The seasonal component showed no abrupt changes, and the seasonal amplitude was approximately 0.05 NDVI.   [95] which describe the water resources state from more to less abundant for: normality (green), pre-alert (yellow), alert (orange), and emergency (red). Table 1 summarizes the date of abrupt changes, their magnitude and direction, and the regression parameters of the piecewise model, which provides information about the NDVI on the break date (intercept) and on the recovery rate (slope). The timing of occurrence and the direction of all four abrupt changes were similar for both treatment and control zones, since their breakpoint dates (95% CI) overlapped. All identified breakpoints occurred immediately after environmental disturbances (i.e., anomalies in water resources in this study). In the first and second abrupt changes an increase in NDVI was observed associated with normal hydric states (HDI) arising after recurrent states of alert or emergency. Conversely, the third and fourth abrupt change occurred along with an emergency or alert state just after the observed NDVI values moved between normality and pre-alert states of water resources. The magnitudes observed at the four breakpoints were always higher for treatment zone than for the control zone (Table 1). Table 1. Timing and NDVI magnitudes of breaks in the trend component of Sa Dragonera (treatment zone) and Mallorca (control zone) time series. The left and right limits of the breakpoint dates indicate a 95% confidence interval of date estimation. 1 = Timing indicated by year and month. Data format: Year (Month). 2: Difference between the NDVI value at the end of a linear model and the intercept of the next one. 3: Linear regression parameters of segments after the breakpoints date, where the intercept represents the starting NDVI, and the slope represents the greening (+) or browning (-) ratio of the vegetation. The parameters of the starting segment are omitted in this Likewise, the recovery rates (Slope in Table 1) were of a similar level in all four modeled linear regressions, highlighting that these rates are not higher in the modeled segments of the treatment zone after the rat eradication campaign date.
Together, the BFAST approaches consistently describe the dynamics of the vegetation and show similar shift dates of NDVI between zones, as well as the direction of the changes, and trends of the piecewise models. These results indicate that the variability of rainfall and water reserves were the main factors controlling the vegetation dynamics following the rat eradication campaign as observed in its previous data and in the control zone data. Furthermore, in terms of the magnitude of the changes and the recovery rates observed, there was no relationship between the elimination of rats and the processes of resistance or recovery of the studied vegetation of the islet.

Discussion
The presented results indicate that the 2011 rodent eradication campaign in Sa Dragonera Islet (treatment zone) had a negligible effect on the primary production dynamics detected at the Landsat resolution scale. The observed changes in NDVI before and after the eradication campaign appear to be driven by water-related environmental events, which could be deemed characteristic of a Mediterranean system [43,53,55,97]. Positive changes detected are followed by segments with a negative or almost zero recovery rate. These patterns coincide with the transition of HDI from alert or pre-alert prolonged states to states of normality. These vegetation dynamics arise in regions with prolonged periods of water limitation interrupted by wet episodes [54] that may trigger a fast germination of short-lived plants and subsequent increase in NDVI followed by a gradual decrease [50]. Other papers such as [59] have used the visual comparison of the rain patterns with dates of abrupt changes of NDVI to evaluate the BFAST algorithm facing known flood events. In this study, we do not validate the dates of the abrupt changes observed due to hydrological factors, since we do not have enough temporal precision in the dates that categorize the HDI. However, the abrupt and gradual change patterns in the NDVI time series of both treatments is used to reject the hypothesis that deratization drove changes in NDVI.
The comparison between the timing of eradication campaign with the timing of the first break following, indicated that BFAST was unable to detect abrupt changes in vegetation response caused by rat elimination. The gradual NDVI increase shown in the fourth segment in the treatment time series (Figure 7a) could reflect the recovery of the ecosystem. However, the BFAST algorithm indicated a non-significant linear relationship, which may have been a relic of noise rather than the signal [59], resulting in the detection of this abrupt change. It is not believed that this represents the response of the vegetation to the demise of rats. The choice of h parameter influences the decomposition of the time series, affecting the change detection process, which is based on signal-to-noise ratio. It is possible that the low h parameters caused the linear models to fit the noise rather than the signal, resulting in the detection of this abrupt change but being non-significant verifies our interpretation. We considered the optimal value of h = 1/7 following the recommendations of [59]. Values of h lower than 1/8 generated additional changes with non-significant linear fits (for example, it identified negative abrupt breaks in 2011), and values greater than 1/6 generated longer segments, which limited the level of detection detail (here it removed changes detected at the end of 2001). No seasonal changes were detected, however Seasonal-Trend decomposition (STL) showed that a seasonal component was present in the time series, therefore it was convenient to remove seasonal effects before analyzing the long-term trends.
The BFAST method dates and quantifies the changes in vegetation from global to local scale, which can be geolocated pixel by pixel at the satellite spatial resolution [43,55]. Often studies using BFAST validate results by prior knowledge of the study area or by subsequent field visits. Our study, however, has been carried out in an entirely remote manner where the NDVI refer to averaged pixels from the sample plots selected under strict spatial criteria, rather than to geolocated pixels. Thus, changes detected characterize the dynamics on the targeted area over time. Averaged pixels allow comparison of the time series of the study area with that of a control zone, where sampling plots that follow similar dynamics could be selected. By doing this the study overcomes the lack of information available about the previous state of the rodenticide campaign.
The presented results do, however, demonstrate that using BFAST is a suitable approach to detect abrupt and gradual changes in NDVI, for example due to the effects of rodenticide campaigns which are known to lead to recovery of vegetation productivity. Within the treatment zone the BFAST model dates and quantifies the NDVI abnormal behavior, considering seasonal variations In this study, we cannot conclusively attribute NDVI recovery, detectable on a 30m spatial scale across the treatment zone, to the eradication campaign, particularly as meteorological drivers in such settings are evident as the dominant influence.
The recovery of an island ecosystem invaded by rats must necessarily begin with their eradication. However, the time and rate of recovery depends on multiple factors such as the frequency, severity and duration of the impacts. There are some studies that have reported the recovery of the vegetation cover in deratized islands through field work, but in an extension and magnitude appreciable by the Landsat sensors. For example, eight years after the deratization of a small island in the Indian Ocean, they reported, through field work, an increase in vegetation cover (from 30% to 70%) validated through aerial images. However, it was due only to the regeneration of the herbaceous stratus [29]. In a study to assess rodent eradication in the Montebello archipelago, Western Australia, vegetation density was monitored for 26 years with Landsat images [34], and the authors reported a trend of vegetation recovery on the islands just two years after rat eradication. Interestingly, and converse to the findings of our study, they observed that the positive relationship typically exhibited between precipitation and vegetation cover was nullified in the rat-infested islands [34].
That said, in the case of Sa Dragonera Islet, some signs of recovery have been described such as the presence of new plant species, the greater presence of other plants that were previously rare, or the greater abundance of arthropods [78].
The original hypothesis suggests the greening of vegetation, especially considering that plants can keep sprouts (and hence NDVI increment) after rat browsing demise. However, our results reflect the lack of change to primary production.
It appears likely that the lack of rat herbivory pressure that could lead to higher plant recruitment and increase in existing plant biomass, did not translate into an abrupt change in the NDVI time series; and that the NDVI values recorded in the Mediterranean ecosystem, before and after the eradication campaign, are principally driven by waterrelated environmental events. The remote sensing record would therefore offer no signs of long-term recovery.
The vegetation of Sa Dragonera Islet and Mallorca Island exhibit a long history of permanent pressures by natural and domestic animals plus recurrent natural and humaninduced fires [79,80]. The herbivorous mammals such as Myotragus balearicus lived in Baleares archipelago till its extinction shortly before the introduction of goats 4000 years ago [98]. Feral goats are currently common in Mallorca Island and they were also in Sa Dragonera Islet until its eradication in 1975 [69,70]. Thus, it could be suggested that the prolonged coexistence of the dominant plant taxa with herbivores in Baleares archipelago, might explain the resistance of vegetation dynamics to herbivory by rats in Dragonera Islet.
Ecosystem recovery is complex and involves cascading processes, and, when it has been severely impacted, recovery could take thousands of years [99,100]. Thus, the absence of remotely detectable NDVI recovery in the treatment zone may be due to the action of disturbance within an environment with inherently stressful conditions (e.g., xeric environment). The presence of black rats could go back several centuries [17] and contribute to disturbance through extremely high population densities. That, combined with historical anthropogenic activity make the islet ecosystem particularly susceptible to land degradation [101]. Therefore, due to an intense and prolonged disturbance the ecosystem compositions could have drastically changed, resulting in a limited recovery following removal of an early driver of degradation in the location [17,30,102,103].
Alternatively, the null response of the vegetation could be due to the fact that the pressure exerted by rats on the plants is not sufficiently high to impact the primary productivity of the island and therefore the recorded NDVI. It could be argued that the ecosystem tolerates the rat herbivore in an already naturally stressed environment [67,104], with most of the dominant plants in the landscape also being sprouting species adapted to browsing and associated disturbance [65][66][67]104]. To assume or discard this option, it would be necessary to have more detailed information on the situation prior to the arrival of the rats on the islet. Likewise, more detailed data on the influence of rats on the richness of plant species, growth rates and periods, seed production, and the consequent population recruitment on the islet would prove vital to deeper functional understanding.
Although records of the ecological past of the islet are sparsely available and limited, the evolution of the ecosystem following elimination of the rats can be monitored in this manner to determine long term trends. The effects of rats on the structure of vegetation are sufficiently reported in the literature, however monitoring eradication outcomes is sporadic and limited [32,35,40], and pre-management data is usually lacking [38]. Hence, this satellite-driven approach represents a consistent methodology that is critical to understanding the ongoing and evolving dynamics of island ecosystems following rat eradication.

Conclusions
This study is the first of its kind to examine the effectiveness of remote sensing and BFAST timeseries analysis in island ecosystems post black rat eradication. Given the general lack of reports assessing the effect of rat eradication on island ecosystems, a satellite approach to assess changes in vegetation productivity over time is proposed. Addressing this generalized information gap, the BFAST model is found to be an appropriate tool to monitor NDVI time series due to its sensitivity to reflect the long-term and short-term changes in vegetation growth at local scales accurately. To interpret the changes found, a 21year time series of NDVI data extracted from satellite images was analyzed and compared with hydrological data, which was found to be the likely driver of NDVI changes in water limited environments as is the case for our study area. These climatological influences are deemed to mask the detectable effects of rat eradication in remotely sensed NDVI data.
The BFAST results reveal that the primary production in the treatment area is sensitive to water cycles but not immediately to cessation of rat activity. In the short term there was no abrupt increase in vegetation primary productivity as a response to the cessation of rat browsing with this significant ecological event being overshadowed by the state of water emergency that occurred close to the time of the eradication campaign. Nor were changes observed in the long term, since facing the aforementioned environmental stress events, the recovery rates in the treatment zone did not overcome those found in the control zone, rather, they continued to show similar trends. Changes may have only occurred at a small spatial scale, as evidenced by the reported field observations of the anecdotal studies in the area. Thus, given the spatial resolution of the core satellite data these changes had little impact on the averaged NDVI produced over the islet. Unlike other studies, changes were not widespread over the island's surface, so it is not possible to detect them at the scale at which Landsat sensors record. We suggest that the non-response of NDVI to rat eradication on Dragonera Islet is due to a high resistance of vegetation to rat predation achieved during a long period of colonization in an already naturally stressed environment.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to the data size.

Conflicts of Interest:
The authors declare no conflict of interest.