Analysis of Site-dependent Pinus halepensis Mill. Defoliation Caused by ‘Candidatus Phytoplasma pini’ through Shape Selection in Landsat Time Series

Abstract

High levels of ‘Candidatus Phytoplasma pini’ have produced extensive forest mortality on Pinus halepensis Mill forests in eastern Spain. This has led to the widespread levels of forest mortality. We used archival Landsat imagery and shapes algorithm implemented in the Google Earth Engine to explore the potential of the LandTrendr algorithm and its outputs, together with field observations, to analyze and predict the health status in P. halepensis stands affected by ‘Candidatus Phytoplasma pini’ in Andalusia (south-eastern Spain). We found that the Landsat time series algorithm (LandTrendr) has captured both long- and short-duration trends and changes in spectral reflectance related to phytoplasma disturbance in the Aleppo pine forest stands investigated. The normalized burn ratio (NBR) trends were positively associated with environmental variables: Annual precipitation, mean temperature, soil depth, percent base saturation and aspect. Environmental variables were tested for their contributions to the mapping of changes in Aleppo pine cover in the study area, as an empirical modeling approach to disturbance mapping in forests of south-eastern Spain. The methodology outlined in this paper has produced valuable results that indicate new possibilities for the use in forest management of remote-sensing technologies based on spectral trajectories associated with pest-diseases defoliation. Given the likely increase in pest risks in the forests of southern Europe, accurate assessment and map of pest outbreaks on forests will become increasingly important, both for research and for practical applications in forest management.

1. Introduction

Worldwide, forest decline and tree mortality events due to climate change have increased since the mid-20th century [1], affecting thousands of hectares in Europe. In the south-western Mediterranean Basin, many coniferous forests have experienced recent large-scale tree mortality events, affecting Scots pine (Pinus sylvestris L.) [2]; black pine (Pinus nigra Arnold.) [3]; mixed conifer stands (Abies pinsapo Boiss.—Pinus spp.) [4] and maritime pine (Pinus pinaster Aiton.) [4], among others [5]. The causes of forest mortality and the conditions that give rise to outbreaks of widespread mortality are complex and hard to predict [6], but emergent forest pests cause regional-scale tree mortality [7].

Many recent studies have emphasized that pest-caused tree mortality related to climate change may play an important role in the Mediterranean Basin [1,2,8]. For instance, [8] evaluated crown defoliation patterns and trends during a 20-year period in a large area of Europe, concluding that defoliation increased only in southern Europe. Similarly, Pinus forest mortality linked to climatic factors has been documented in southern Spain [2,4,9] although more-proximal exogenous factors, such as forest pests, and endogenous factors, such as stand density and environmental setting, are also important [10]. Recently, ‘Candidatus Phytoplasma pini’ has been nominated as a potential severe pest of pine species, including Pinus halepensis Mill. The symptoms are characterized by shoot proliferations and needle aberration [11]. The phytoplasma detected is widely dispersed in many different Aleppo pine forests in southern Spain, with important environmental and economic impacts.

The monitoring and description of these mortality processes are complex tasks since there are often multiple stressors causing a certain symptom [12], so the key problem of understanding and mapping such pest interactions needs to be overcome. The monitoring of this decline and mortality has occurred largely through the use of multi-temporal and multi-spectral remote sensing techniques, which are powerful tools to detect and study vegetation changes in an inexpensive, fast and periodical way [13,14]. Previous Landsat temporal data have been employed to map forest pest disturbances at a variety of spatial and temporal scales [15,16]. Landsat is the only remote sensing mission that offers more than 30 years of continuous spectral information at medium spatial resolution (pixels of 30–120 m); this enables the effective study of vegetation disturbance and changes in spectral values related to pest impacts [17]. In spectral terms, vegetation, when healthy, shows the highest absorption and the lowest reflectance in the visible region, and the other way round in the near infrared region (NIR) [18]. Differences in the spectral response due to vegetation health status have resulted in the use of vegetation indices to monitor stress. Indeed, various vegetation indices, such as Normalized Difference Vegetation Index (NDVI), have been widely used to monitor changes in vegetation. Another index, the normalized burn ratio, which contrasts the Near Infrared (NIR) and Panchromatic and Shortwave Infrared (SWIR) bands, has been widely used to track the effects of pests on vegetation [15,19,20,21]. More recently, Landsat-based algorithms have emerged for the detection of the changes in the short-term anomalies of spectral trajectories related to the effects of defoliation caused by pests. These include the vegetation change tracker (VCT), vegetation continuous fields (VCF), LandTrendr, continuous change detection and classification (CCDC), multi-index integrated change algorithm (MIICA), a Fourier regression algorithm, a gradual ecosystem change algorithm and methods based on non-parametric statistical analysis [22]. Therefore, the relationship between the spectral trajectories of Landsat image time series and the temporal dynamics of pest disturbance can be applied to detect subtle and abrupt inter-annual forest changes [13,23]. The combination of free and temporally dense data from Landsat and temporally based algorithms opens the door to the mapping of defoliation impacts in a consistent way across large areas of forest that are affected, around the world [24,25]. In addition, maps of the defoliation due to pests are required to understand the interaction between the potential pest damage and key biophysical drivers (e.g., the physical structure of the stand, and climatic and biotic disturbances) [14].

In P. halepensis forests of south-eastern Spain, ‘Candidatus Phytoplasma pini’ outbreaks can induce relatively rapid tree defoliation and mortality associated with typical symptoms including proliferation of axillary buds (resulting in a witches’ broom growth), abnormal internodes elongation and generalized stunting and eventual total tree crown mortality [26]. This pest-caused tree mortality necessitates alterations in forest management and eventually increases the risks of forest fires and severe erosion. However, there is a lack of studies involving the mapping of pest-related defoliation in Mediterranean forests based on Landsat time series change detection algorithms [27] and ground-based measurements.

In this paper, we explore the potential of the LandTrendr algorithm and its outputs, together with field observations, to analyze and predict the health status in P. halepensis stands affected by ‘Candidatus Phytoplasma pini’ in Andalusia (south-eastern Spain). In relation to this, we had three specific aims: (i) To characterize spectral trajectories in Landsat time series associated with Candidatus Phytoplasma pini defoliation disturbances of varying duration and severity; (ii) to relate changes in the normalized burn ratio (NBR) index spectral trajectories of the last segment (∆NBRls) to environmental variables, to identify the major drivers of phytoplasma-caused defoliation and tree mortality through a the k-nearest neighbors algorithm (k-NN) and (iii) to use a model to predict the vulnerability of the Aleppo pine stands in our area of study. In this way, we were able to study and map how environmental factors affect the phytoplasma-related defoliation and mortality rates of P. halepensis, which determine the potential management alternatives in a changing environment.

2. Materials and Methods 2.1. Insect Outbreak Area

We focused this study on Pinus halepensis forests in the south-eastern Iberian Peninsula that have experienced contrasting phytoplasma (CPP) activity, according to the Forest Health Andalusia Network (Red SEDA, Figure S1 Supplementary Material), beginning in 2013 and increasing with time, following the sporadic and locally intense drought-induced decline processes in the early 2000s. Symptoms associated with ‘Candidatus Phytoplasma pini’ on leaves and twigs of P. halepensis were observed and a positive diagnosis of phytoplasma presence was verified. The study area was located in southeastern Spain (Guadix-Baza, Granada, hereafter SGB, 37°21′ N 3°5′ W, elevation 750 m to 2270 m, Figure S1 Supplementary Material). The climate is Mediterranean semi-arid, with estimated annual precipitation for the period 1990–2016 ranging between 354 mm and 494 mm, concentrated mainly in winter and to a lesser extent in spring–fall. The average annual temperature is between 10 and 15 °C, highlighting the marked contrast between the winter months (monthly average close to 5 °C) and the summer months (monthly average around 25 °C). The average annual potential evapotranspiration (PET) is more than twice the average annual rainfall, the water balance having negative values for 7–8 months of the year. The soils have a low water holding capacity, with the most-abundant soil types being entisols and inceptisols on limestones. The landscape is characterized by the presence of badlands originated by the erosive action of water on materials of certain plasticity such as clays, sands and silts. P. halepensis plantations in this area are characterized by high stand density because of intensive forest plantations in the 1960s and a subsequent lack of silvicultural practices (Table S1 Supplementary Material).

2.2. Field Data

In November 2016, we placed 72 survey plots (60 m × 60 m, each comprising the extent of four Landsat pixels) across a range of phytoplasma defoliation levels, but with similar forest structure conditions. Forest mortality patterns were quantified across the study area using the inventory form from the forest management plan (2009, Andalusia Forest Service; Figure S1, Supplementary Material). We excluded all non-forested areas, areas that had received silvicultural treatments, areas dominated by shrub cover, areas dominated by other pine species (Pinus pinaster Aiton., P. nigra Arnold. and P. sylvestris L.), riparian vegetation and areas within 60 m of major roads. We also verified, using an orthophoto of the area, that the pixels were completely inside the P. halepensis stands, and that other typologies (such as paths or roads) were not included in the pixels clusters. In each plot, we measured the defoliation level—removal of photosynthetically active leaves—according to [28] the basal area (m² ha⁻¹) and mortality, as the ratio between the living and dead trees inside the plot. All plots were located in an area having a previous positive diagnosis of the presence of ‘Candidatus Phytoplasma pini’. These plots were classified in three damage levels: ‘‘High” (defoliation-mortality covered 60–100% of the polygon), ‘‘Medium” (defoliation-mortality covered 26–59% of the polygon) and ‘‘Low” (defoliation covered <25% of the polygon, with no evidence of mortality; Figure S2 Supplementary Information). Besides the field plots, we set up 600 plots over the area cover by the P. halepensis stands, delimited by the distribution of Aleppo pine in SGB. The plots were distributed randomly by using the random points inside the polygons tool in QGIS 2.18 [29], with a minimum distance between points of 100 m. After the cleaning of plots placed outside forested areas or near to roads, the total number of plots was 575.

2.3. Image Pre-Processing and Vegetation Index Calculation

Our study used several data sets and required the development of remote sensing indices and data analysis procedures. Therefore, a flow chart outlining the steps and relationships of each process is provided in Figure 1.

The processing of the Landsat data was carried out using the Google Earth engine’s LandTrendr scripts (LT-GEE ) developed by [27,30]. First, we built the stacks of Landsat Thematic Mapper TM, enhanced thematic mapper + (ETM+) and operational land imager (OLI) imagery for the period 1986 to 2016, using USGS Landsat Surface Reflectance Tier 1 data sets.

The LT-GEE application programming interface (API) includes a function to create the annual composites required for the LandTrendr algorithm, buildSRcollection. This function builds an annual cloud and cloud shadow masked, created by CFMASK (C code based on the Function of Mask Fmask algorithm) [31], and a composite generated by a medoid selection process [32]. OLI image bands 2, 3, 4, 5, 6 and 7 were a subset and were transformed to the spectral properties of ETM+ bands 1, 2, 3, 4, 5 and 7, respectively, using slopes and intercepts from reduced major axis regressions reported in [33]. To minimize the presence of clouds and the effects of vegetation phenology, sun angle differences and other factors, only images gathered from 15 June to 30 September were selected.

We then calculated for each pixel the normalized burn ratio (NBR: [NIR-SWIR2]/[SWIR2+NIR]; [20,34], to build the pixel-level time series through the transformSRcollection function.

2.4. LandTrendr Outputs and Health Status of the Plots

Image NBR stacks were processed with the LT-GEE segmentation algorithm in the GEE user interface. The parameters used in the algorithm are shown in Table S2 (Supplementary Material). The output is delivered as a GEE image with three bands: One containing annual segmentation information; another including the spectrally fitted annual series (FTV) values and the last one with the Root Mean Square Error (RMSE) of the segmentation fit.

We extracted the segmentation information by using the getSegmentData function. The result is an image array with dimensions: 8 (rows) × nSegments (cols). Each row describes an attribute of the segments identified by LT per pixel time series and this can be used to consult the meaning of each row, see https://emapr.github.io/LT-GEE/api.html#getsegmentdata. Each column represents a segment in the time series per pixel, ordered from the earliest to the latest year. Disturbance segments have been defined as those segments experiencing a decline in NBR over time [27]. We selected the spectral delta corresponding to the last segment value of the fitted series (hereinafter, ∆NBR_ls) to evaluate the health status of the stands. Finally, we established three categories of health status according to the ∆NBR_ls values obtained in the sample field plots and their defoliation level. An analysis of variance (ANOVA) procedure was used to determine if the defoliation levels measured in the field plots differed significantly among the health categories.

2.5. Environmental Variables Related to ∆NBRls

The last step was to assess the influence of the site conditions on the ∆NBR_ls values, in order to predict the health status of the pine stands through a k-nearest neighbor algorithm. Therefore, we obtained bioclimatic variables, at a spatial resolution of 200 m × 200 m, compiled in the Andalusia Environmental Data Network (REDIAM, available from https://descargasrediam.cica.es/repo/s/RUR), for the study area. The dataset (40 variables) contains four categories of variables: 16 climatic variables, 12 edaphic variables, 12 topographical variables and four tree cover variables (Table S3, Supplementary Material). Climate data were obtained from the meteorological station of Baza (Regional Climate Centre, 1960–2016, 37°33′52″ N, 02°46′03″ W, 814 m). The topographical variables were generated from a digital elevation model (5 m spatial resolution) obtained from the REDIAM, in QGIS version 2.18 [29]. These variables were chosen based on their eco-physiological relevance to the tree drought response [4,35].

Prior to development and implementation of the k-NN algorithm, the following steps were taken to prepare the proposed environmental variables for analysis. First, a variance inflation factor (VIF) analysis was carried out among the predictor variables to determine the presence of multicollinearity [36]. Variables with a VIF >10 were removed from the posterior analysis. The collinearity analysis was performed in R 3.5.1 [37] using the R package usdm [38]. Second, we selected the most-important variables that were not highly correlated for the imputation modeling, using the varSelection function from yaImpute [39]. This function computes the generalized root mean square distance (GRMSD) as variables are added to a k-NN imputation algorithm. By adding variables, the varSelection function keeps variables that strengthen the imputation and excludes variables that weaken the imputation the most [39]. Finally, we performed the k-NN imputation (k = 3) to characterize the relationship between the environmental predictor (selected environmental variables) and the response (∆NBR_ls) variables. We used random forest (RF) [40].

In order to validate the estimations of the models, the input data were separated into training and evaluation sets, covering 70% and 30% of the input data respectively. Later, the variable selection process was executed by choosing the combination of variables that minimized the GRMSD when variables were added or deleted one at a time [39]. After the best model was chosen, internal validation to assess the model accuracy was implemented, including a Q-value overfitting test. External validation and cross-validation were then run. The RMSE and BIAS (1 and 2) were calculated for each process.

RMSE=_∑i=1N ^{(_yi^−_yi)2}n

BIAS=_∑i=1N(_yi^−_yi)n

where y_i is the observed value, ŷ_i is the estimated value and n is the number of observations.

2.6. Map of the Current and Potential Risk of ‘Candidatus Phytoplasma pini’

Finally, we generated a map of predicted ∆NBR_ls attribution classes at a 200 m × 200 m resolution for the P. halepensis plantations in the study area. The attribution classes included three levels of potential risk of ‘Candidatus Phytoplasma pini’, according to the values of the predictive environmental variables associated with the current defoliation-mortality decline. Strictly, the purpose of this map was not to be a final mapping product but rather to emphasize both disturbances probably caused by phytoplasma and the potential of Landsat time series data to map pest impacts as a key tool for forest management in high-risk pest areas. For this, we used the predict function from the raster package in R [41].

3. Results 3.1. ∆NBR_ls of Phytoplasma-Affected Pine Stands

In general, phytoplasma effects were evident in the Landsat time series as combinations of both short- and long-duration spectral changes. The spectral trajectories of Aleppo pine NBR produced a particular pattern of shapes output for each of the three damage levels (Figure 2). In fact, the NBR values exhibited long-duration increases (1986–2013) for the three damage levels, corresponding to relatively slow growth processes at an annual time scale. However, stands with moderate and severe levels of damage showed outstanding decreases in ∆NBR_ls trend between 2013 and 2016. The ANOVA showed that the defoliation levels were significantly different among the health categories, according to the ∆NBRls values (F = 6.191; p < 0.001).

3.2. Environmental Predictors of ‘Candidatus Phytoplasma pini’ Defoliation

Investigation of how specific ∆NBR_ls damage estimates respond as a function of the environmental conditions may help to predict potential impacts of ‘Candidatus Phytoplasma pini’. A total of 40 environmental variables were used to perform RF regression, following a sparsity test of the predictor matrix. The ∆NBR_ls damage was related to five environmental variables (Figure 3, Table 1), given their importance in driving ‘Candidatus Phytoplasma pini’ defoliation-mortality. Across the study landscapes, ∆NBR_ls damage was positively correlated (R² = 0.65, RMSE = 0.014%) with variables related to the thermal and precipitation regimes; for instance, the annual precipitation and average mean temperature showed the highest values of importance. The predictor variables related to soil (soil depth and percent base saturation) had also a high-important relationship with ∆NBR_ls damage levels (Figure 3). Figure 4 shows the scatter plots for correlations contrasting the observed versus the estimated values for all ∆NBR_ls damage predictions.

3.3. Disturbance Attribution Maps

Finally, a map of the predicted ‘Candidatus Phytoplasma pini’ defoliation of the P. halepensis forests in the study area was generated, based on the use of the RF model spatial information on the levels of damage to distinguish predicted ∆NBR_ls values as a risk map of ‘Candidatus Phytoplasma pini’ disturbance. Figure 5 illustrates the patterns of the risk of disturbance, showing the potential impact of pest disturbance throughout the P. halepensis forests in the study area. The majority of the area (71.8%) was designated as middle-risk, another 16.3% fell within the low-risk category and 11.8% was within high-risk areas (Figure 6). It should be noted that in the study site 97% of the surface area has been affected already by ‘Candidatus Phytoplasma pini’. In other parts of the study area, dominated by P. halepensis, severe disturbances could be causing large areas of mortality.

4. Discussion

The causes of forest mortality are inherently complex and difficult to predict [42]. In our study, Aleppo pine defoliation and mortality has been related to ‘Candidatus Phytoplasma pini’, generally associated with extensive defoliation processes leading to widespread forest mortality events. Therefore, this paper documents the use of the LandTrendr algorithm and Landsat multi-temporal analysis: (1) To study the feasibility of using the variability in the temporal spectral trajectories of the NBR index to map the current location of phytoplasma defoliation in Aleppo pine forests, and (2) to map the potential risk according to three levels of defoliation as a function of environmental variables. As part of this effort, the ∆NBR_ls shapes trends were analyzed using the Google Earth engine, to provide forest pest-damage mapping as a key tool to develop forest management alternatives. With our Landsat-based approach, the current and potential damage due to phytoplasma were mapped, and their application in adaptive forest management of pests and in detection alarms involving satellite tools could be highlighted.

4.1. Temporal Trends of Remotely Sensed Data

During recent decades, Aleppo pine forests in south-eastern Spain have been prone to physiological changes and dieback processes, which could undergo expansion and contraction processes over time, due to both the pests action, and its combination with climatic changes, particularly severe droughts [43]. In the study area, ‘Candidatus Phytoplasma pini’ has been identified as the major forest pest on this species by excluding other potential causes of biotic or abiotic stresses. In general, Pinus stands affected by phytoplasma exhibit specific foliar symptoms—such as an increased proportion of defoliated trees and lower photosynthetic efficiency—and tree mortality before a more abrupt collapse [44]. Our field observations agreed with other studies that have found severe defoliation and related mortality of pine-dominated forests in response to this pest in European temperate zones [45]. Such changes in the leaves and crowns suggest that remote sensing could be useful for the development of disease-forecasting methods [18].

Among all the available vegetation indexes, NBR has been used frequently to monitor defoliation [21]. In this study, temporal trends of the Landsat series of the NBR index showed differences in values according to the phytoplasma defoliation levels, with different diagnostic disturbance trajectories. Initially, the NBR index showed a progressive long-duration increase in all stands, corresponding to relatively slow growth and non-severe disturbance processes until 2013. Then, the observed differences in the Landsat time series, manifested as short-duration NBR trajectories (variable combinations of both long and short duration), reflect the defoliation and mortality effects of the three types of phytoplasma damage. Indeed, we observed—between 2013 and 2017—a sharp decrease in the NBR values for the medium and high levels of defoliation and mortality, coincident with the period of highest mortality in these forests (SEDA-Network data unpublished). Phytoplasma induces very rapid crown changes and its pixel-scale impacts, expressed as the NBR index, dropped because of changes in chlorophyll content and defoliation in the stands over two or three years. We suggest that the short-term and sharp changes in the tendency of the NBR index to decline are a diagnostic spectral trajectory for phytoplasma in such P. halepensis forests, where part of the canopy is defoliated, resulting in relatively slight re-greening or death. This agrees with some studies that described the effect of defoliators on the spectral response patterns in other pine species [46]. Meigs et al. [7] analyzed the defoliation of pine forests in the US Pacific Northwest Region, using different vegetation indexes, and they observed that the trend of the NBR magnitude distinguished insect disturbance, since it was sensitive to the increasing predominance of conditions that led to defoliation. The average values of NBR of symptomatic stands can be related to symptoms severity and interannual changes. In contrast, a continuous increase in the NBR value was observed in asymptomatic stands, therefore, NBR can be used to discriminate the “healthy” and the “affected” forests.

Analysis of the temporal trends of NBR images could allow prompt identification of the early stages of the symptoms, even before they can be visually detected. However, the NBR measurements have clear limitations, because they provide values capable of identifying alterations in photosynthetic efficiency, but do not discriminate the causes of biotic (pest and diseases) or abiotic changes (e.g., silvicultural treatments) or stress (e.g., droughts) afflicting plants. Our work was possible because the field plots were selected to avoid other causes of changes in canopy defoliation, confirming that alterations in NBR were exclusively linked to symptoms of phytoplasma. Early detection of pests, before their symptoms appear, could be integrated into the planning of silvicultural treatments for forests. However, no single model of temporal change can capture the full range of pest-induced changes across heterogeneous landscapes, so, remote-sensing approaches still require the integration of both spectral trends, to detect the full variety of short- and long-duration changes, and field observations to be reliable predictors of forest pest damage [27].

4.2. Relationship between NBR and Environmental Variables

Across the P. halepensis forests studied, the LandTrendr disturbance magnitude (expressed as ∆NBR_ls) was positively correlated with some environmental variables according to a kNN algorithm (R² = 0.62, p < 0.0001). Climate (annual precipitation, average mean temperature and average net primary production) and soil (percent base saturation and soil depth) were identified as the most important environmental factors related to phytoplasma-induced defoliation in this study. P. halepensis forests located on shallow soils, with lower rainfall, facing south and with thermal (higher temperatures) and edaphic limitations (more basic soils) were shown to suffer greater impacts of phytoplasma. In such locations drought-related stress may be greater despite the strong drought-tolerance characteristics of Aleppo pine—which likely exhibits stomatal control of water loss and cavitation resistance, as has been shown in xeric environments [47]. Precipitation and thermal parameters are known to have a crucial influence on the vegetative response of Aleppo pine forests [48]. In 2013 there was very little rainfall in the study area, and it is probable that leaf metabolism was affected and predisposed to pest damage. These physiographic factors seem to be synergic with other, site-specific environmental conditions—such as soil depth or fertility (base saturation)—that modulate the drought response of Aleppo pine. By contrast, locations with water stress compensation experienced fewer forest defoliation and mortality events and these were both more episodic and smaller in magnitude. Our results show that the importance of phytoplasma-related defoliation-mortality has been greatest for drier forest locations and during drought periods, a phenomenon modulated by the effect of soil conditions closely linked to climatic variability. These results accord with numerous studies of pest-forest interactions that found more severe damage in environments with high abiotic stress and competition (e.g., [49]) and are consistent with our understanding of the biological response of phytoplasma. The importance of these environmental variables could be modified by the forest structure, in particular tree density, which was not considered in this study. Denser stands may facilitate the spread of pathogens [50] or may indicate previously high densities used in plantations on microsites that were or currently are unfavorable. Low-density and high-quality sites at wetter microsites can ameliorate the combination of pests and harsh environmental extremes, but the negative effect of drought periods and dense stands may play an important accelerating role, favoring defoliation and mortality.

These results seem to reinforce the idea that an appropriate use of field data together with LandTrendr could provide a more objective, consistent and quantitative measure of the magnitude of phytoplasma defoliation-mortality once the damage has already appeared. Nevertheless, the use of LandTrendr could allow us to anticipate the damage associated with phytoplasma before the appearance of visual symptoms such as chlorosis or defoliation, because the field plots cover the full gradient of damage across highly heterogeneous stands. Future research could investigate the effect of density on phytoplasma stand interactions. Additionally, the use of better spatial resolution environmental information (e.g., more appropriate to the spatial resolution of the sensors used) could improve the identification and interpretation of the environmental drivers involved in the woodland mortality processes.

4.3. Landsat as a Source of Phytoplasma Risk Maps and Management Implications

Our previous results, obtained at a medium spatial resolution in Aleppo pine forests, suggest that ΔNBR trends derived from Landsat data are suitable for the mapping of phytoplasma-provoked forest defoliation at the stand and landscape scales. The current and potential risk maps describe a patchy spatiotemporal coverage, including the damage levels classification variability in order to estimate the magnitude of cumulative defoliation-mortality due to phytoplasma, relative to ground-based estimates. As a relative measure of pest effects on P. halepensis stands, these maps captured the cumulative defoliation-mortality spatial variation, showing geographic patterns of damage levels and spatial magnitude and thus providing invaluable identification of specific pest impacts. Due to the cost associated with extensive field data collection, the use of Landsat data, based on a quantitative relationship between field and remote estimates, is particularly encouraging. Current and potential damage map—together with larger field datasets covering multiple time intervals and a full gradient of vegetation-pest conditions, such as the regional Forest Health Network of Andalusia (SEDA-Network)—will provide valuable spatial–temporal coverage throughout the Mediterranean P. halepensis forests, as well as other natural and planted conifer forests. This more detailed information could help the development of management alternatives for vulnerable forests to specific insect agents as well as the determination of pre-visual detectability thresholds.

However, since pest disturbance processes are dynamic and continue to evolve, it should be noted that the discrimination between levels of damage is complex and not always significant between closer levels of damage (e.g., medium and high levels of damage). Some of the limitations are related to the fact that the spectral response of single polygons along the temporal series can be influenced by many factors (including sun angle, phenological variation, the scale of observation and productivity) and by the internal heterogeneity of stands due to heterogeneous defoliation-mortality at the tree scale. Despite these limitations, the use of Landsat-based approaches in the mapping of phytoplasma effects, as a predictive tool, could open new perspectives in the knowledge and control of this important pine disease, as part of the increasing development of precision silvicultural approaches in forest management [27].

5. Conclusions

To the best of our knowledge, no previously published study has reported a direct comparison of field data and Landsat temporal data for tree defoliation-mortality caused by phytoplasma. The methodology outlined in this paper has produced valuable results that indicate new possibilities for the use in forest management of remote-sensing technologies based on spectral trajectories associated with pest-diseases defoliation. A Landsat time series algorithm (LandTrendr) has been used to capture both long- and short-duration trends and changes in spectral reflectance (based on the NBR index) for the mapping of three spectral trajectories that could be diagnostic for phytoplasma disturbance in the Aleppo pine forest stands investigated. The NBR trends were positively associated with environmental variables, showing the potential to develop Landsat-based defoliation maps. As part of this effort, the shapes algorithm was implemented in the Google Earth engine, providing shapes parameters and fitted trajectories for use in phytoplasma risk mapping. Environmental variables were tested for their contributions to the mapping of changes in Aleppo pine cover in the study area, as an empirical modeling approach to disturbance mapping in forests of south-eastern Spain. Our results highlight the utility of pest-mapping methods to cover a wide range of site conditions based on relatively short-term anomalies. Given the likely increase in pest risks in the forests of southern Europe, accurate assessment and map of pest outbreaks on forests will become increasingly important, both for research and for practical applications in forest management.

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

Variables	COD	Low	Medium	High
Annual precipitation (mm)	PRC	470.4 (2.8)	362.4 (2.1)	382.4 (2.3)
Average mean temperature (°C)	T_MED	13.8 (0.5)	15.0 (0.4)	14.3 (0.5)
Average net primary production	DF	1381.8 (31.3)	745.2 (19.7)	988.6 (28.7)
Percent base saturation (%)	PBS	98.7 (0.1)	99.5 (0.03)	99.87 (0.01)
Soil depth (cm)	PS	73.7 (1.4)	71.6 (1.3)	68.3 (1.4)
Topographic exposure (°)	TP_EXPO	–3.45 (7.4)	–1.91 (6.0)	5.43 (6.6)

Supplementary Materials

The following are available online at https://www.mdpi.com/2072-4292/11/16/1868/s1, Figure S1. Pinus halepensis Mill forests location in South-Eastern Spain and detail of the sampled stands in this work, Figure S2. Field defoliation percentage of Pinus halepensis Mill in South-Eastern Spain caused by ‘Candidatus Phytoplasma pini’, Table S1: Forestry characteristics of each species from Third National Forest Inventory of Spain, Table S2. LandTrendr parameters used for the segmentation of the NBR Landsat time series, Table S3. Environmental data used to predict the defoliation of Pinus halepensis Mill in South-Eastern Spain caused by ‘Candidatus Phytoplasma pini’.

Author Contributions

Conceptualization, R.M.N.-C. and J.T.-T.; Methodology, J.T.-T. and R.M.N.-C.; Formal Analysis, J.T.-T. and R.M.N.-C.; Investigation, R.M.N.-C. and J.T.-T.; Resources, R.M.N.-C.; Data Curation, J.T.-T.; Writing-Original Draft Preparation, R.M.N.-C. and J.T.-T.; Writing-Review & Editing, R.M.N.-C. and J.T.-T.; Project Administration, R.M.N.-C.

Funding

This research was funded by the Ministerio de Ciencia, Innovación y Universidades (Spain), through the ESPECTRAMED (CGL2017-86161-R), ECOPOTENTIAL (grant agreement No. 641762), and ISO-Pine (UCO-1265298) projects.

Acknowledgments

The authors thank the Andalusia Department of Agriculture and Environment, which provided access to and background information on the field site, in particular to the forest health technicians (Red SEDA), JM Ruiz Navarro, MA Gómez de Dios, and A Muñoz Risueño. Also to mention the help provided by “Plan Propio de la UGR” and “Programa de Unidades de Excelencia”, for the publication of this article. We are very grateful to Rafael Sánchez, Emidio Silveiro, and Ana Moyano, for their valuable assistance during fieldwork and data acquisition and processing (without them, this project would not have been possible), and to David Walker, for the linguistic revision of the manuscript. We also acknowledge the institutional support of the University of Cordoba-Campus de Excelencia CEIA3.

Conflicts of Interest

The authors declare no conflict of interest.

References

1. Allen, C.D.; Macalady, A.K.; Chenchouni, H.; Bachelet, D.; McDowell, N.; Vennetier, M.; Kitzberger, T.; Rigling, A.; Breshears, D.D.; (Ted) Hogg, E.H.; et al. A global overview of drought and heat-induced tree mortality reveals emerging climate change risks for forests. For. Ecol. Manag. 2010, 259, 660–684.

2. Sánchez-Salguero, R.; Navarro-Cerrillo, R.M.; Camarero, J.J.; Fernández-Cancio, Á. Selective drought-induced decline of pine species in southeastern Spain. Clim. Chang. 2012, 113, 767–785.

3. Sánchez-Salguero, R.; Camarero, J.J.; Dobbertin, M.; Fernández-Cancio, Á.; Vilà-Cabrera, A.; Manzanedo, R.D.; Zavala, M.A.; Navarro-Cerrillo, R.M. Contrasting vulnerability and resilience to drought-induced decline of densely planted vs. natural rear-edge Pinus nigra forests. For. Ecol. Manag. 2013, 310, 956–967.

4. Navarro-Cerrillo, R.M.; Rodriguez-Vallejo, C.; Silveiro, E.; Hortal, A.; Palacios-Rodríguez, G.; Duque-Lazo, J.; Camarero, J.J. Cumulative Drought Stress Leads to a Loss of Growth Resilience and Explains Higher Mortality in Planted than in Naturally Regenerated Pinus pinaster Stands. Forests 2018, 9, 358.

5. Peña-Gallardo, M.; Vicente-Serrano, S.M.; Camarero, J.J.; Gazol, A.; Sánchez-Salguero, R.; Domínguez-Castro, F.; El Kenawy, A.; Beguería-Portugés, S.; Gutiérrez, E.; De Luis, M.; et al. Drought Sensitiveness on Forest Growth in Peninsular Spain and the Balearic Islands. Forests 2018, 9, 524.

6. Bentz, B.J.; Régnière, J.; Fettig, C.J.; Hansen, E.M.; Hayes, J.L.; Hicke, J.A.; Kelsey, R.G.; Negrón, J.F.; Seybold, S.J. Climate Change and Bark Beetles of the Western United States and Canada: Direct and Indirect Effects. BioScience 2010, 60, 602–613.

7. Meigs, G.W.; Kennedy, R.E.; Gray, A.N.; Gregory, M.J. Spatiotemporal dynamics of recent mountain pine beetle and western spruce budworm outbreaks across the Pacific Northwest Region, USA. For. Ecol. Manag. 2015, 339, 71–86.

8. Carnicer, J.; Coll, M.; Ninyerola, M.; Pons, X.; Sánchez, G.; Peñuelas, J. Widespread crown condition decline, food web disruption, and amplified tree mortality with increased climate change-type drought. Proc. Natl. Acad. Sci. USA 2011, 108, 1474–1478.

9. Navarro-Cerrillo, R.; Varo, M.A.; Lanjeri, S.; Clemente, R.H. Cartografía de defoliación en los pinares de pino silvestre (Pinus sylvestris L.) y pino salgareño (Pinus nigra Arnold.) en la Sierra de los Filabres. Rev. Ecosistemas 2007, 16, 163–171.

10. Netherer, S.; Schopf, A. Potential effects of climate change on insect herbivores in European forests—General aspects and the pine processionary moth as specific example. For. Ecol. Manag. 2010, 259, 831–838.

11. Schneider, B.; Torres, E.; Martín, M.P.; Schröder, M.; Behnke, H.-D.; Seemüller, E. ‘Candidatus Phytoplasma pini’, a novel taxon from Pinus silvestris and Pinus halepensis. Int. J. Syst. Evol. Microbiol. 2005, 55, 303–307.

12. Trumbore, S.; Brando, P.; Hartmann, H. Forest health and global change. Science 2015, 349, 814–818.

13. Czerwinski, C.J.; King, D.J.; Mitchell, S.W. Mapping forest growth and decline in a temperate mixed forest using temporal trend analysis of Landsat imagery, 1987–2010. Remote Sens. Environ. 2014, 141, 188–200.

14. Hall, R.J.; Castilla, G.; White, J.C.; Cooke, B.J.; Skakun, R.S. Remote sensing of forest pest damage: A review and lessons learned from a Canadian perspective. Can. Entomol. 2016, 148, S296–S356.

15. Wulder, M.A.; Dymond, C.C.; White, J.C.; Leckie, D.G.; Carroll, A.L. Surveying mountain pine beetle damage of forests: A review of remote sensing opportunities. For. Ecol. Manag. 2006, 221, 27–41.

16. Dennison, P.E.; Brunelle, A.R.; Carter, V.A. Assessing canopy mortality during a mountain pine beetle outbreak using GeoEye-1 high spatial resolution satellite data. Remote Sens. Environ. 2010, 114, 2431–2435.

17. Cohen, W.B.; Yang, Z.; Stehman, S.V.; Schroeder, T.A.; Bell, D.M.; Masek, J.G.; Huang, C.; Meigs, G.W. Forest disturbance across the conterminous United States from 1985–2012: The emerging dominance of forest decline. For. Ecol. Manag. 2016, 360, 242–252.

18. Silva, C.R.; Olthoff, A.E.; de la Mata, J.A.D.; Alonso, A.P. Remote monitoring of forest insect defoliation. A review. For. Syst. 2013, 22, 377–391.

19. Hayes, D.J.; Sader, S.A. Comparison of Change-Detection Techniques for Monitoring Tropical Forest Clearing and Vegetation Regrowth in a Time Series. Photogramm. Eng. Remote Sens. 2001, 67, 1067–1075.

20. Key, C.H.; Benson, N.C. Landscape assessment: Remote sensing of severity, the normalized burn ratio and ground measure of severity, the composite burn index. In FIREMON: Fire Effects Monitoring and Inventory System; USDA Forest Service, Rocky Mountain Res. Station: Ogden, UT, USA, 2005.

21. Housman, I.W.; Chastain, R.A.; Finco, M.V. An Evaluation of Forest Health Insect and Disease Survey Data and Satellite-Based Remote Sensing Forest Change Detection Methods: Case Studies in the United States. Remote Sens. 2018, 10, 1184.

22. Moisen, G.G.; Meyer, M.C.; Schroeder, T.A.; Liao, X.; Schleeweis, K.G.; Freeman, E.A.; Toney, C. Shape selection in Landsat time series: A tool for monitoring forest dynamics. Glob. Chang. Biol. 2016, 22, 3518–3528.

23. Vogelmann, J.E.; Xian, G.; Homer, C.; Tolk, B. Monitoring gradual ecosystem change using Landsat time series analyses: Case studies in selected forest and rangeland ecosystems. Remote Sens. Environ. 2012, 122, 92–105.

24. Wulder, M.A.; Masek, J.G.; Cohen, W.B.; Loveland, T.R.; Woodcock, C.E. Opening the archive: How free data has enabled the science and monitoring promise of Landsat. Remote Sens. Environ. 2012, 122, 2–10.

25. Kennedy, R.E.; Andréfouët, S.; Cohen, W.B.; Gómez, C.; Griffiths, P.; Hais, M.; Healey, S.P.; Helmer, E.H.; Hostert, P.; Lyons, M.B.; et al. Bringing an ecological view of change to Landsat-based remote sensing. Front. Ecol. Environ. 2014, 12, 339–346.

26. Dickinson, M.; Tuffen, M.; Hodgetts, J. The Phytoplasmas: An Introduction. In Phytoplasma: Methods and Protocols; Methods in Molecular Biology; Dickinson, M., Hodgetts, J., Eds.; Humana Press: Totowa, NJ, USA, 2013; pp. 1–14. ISBN 978-1-62703-089-2.

27. Kennedy, R.E.; Yang, Z.; Cohen, W.B. Detecting trends in forest disturbance and recovery using yearly Landsat time series: 1. LandTrendr—Temporal segmentation algorithms. Remote Sens. Environ. 2010, 114, 2897–2910.

28. Ferretti, M. Mediterranean Forest Trees a Guide for Crown Assessment; Commission of the European Communities: Brussels, Belgium, 1994.

29. QGIS Development Team. QGIS Geographic Information System; Open Source Geospatial Foundation Project; 2018.

30. Kennedy, R.; Yang, Z.; Gorelick, N.; Braaten, J.; Cavalcante, L.; Cohen, W.; Healey, S.; Kennedy, R.E.; Yang, Z.; Gorelick, N.; et al. Implementation of the LandTrendr Algorithm on Google Earth Engine. Remote Sens. 2018, 10, 691.

31. Zhu, Z.; Woodcock, C.E. Object-based cloud and cloud shadow detection in Landsat imagery. Remote Sens. Environ. 2012, 118, 83–94.

32. Flood, N. Seasonal Composite Landsat TM/ETM+ Images Using the Medoid (a Multi-Dimensional Median). Remote Sens. 2013, 5, 6481–6500.

33. Roy, D.P.; Kovalskyy, V.; Zhang, H.K.; Vermote, E.F.; Yan, L.; Kumar, S.S.; Egorov, A. Characterization of Landsat-7 to Landsat-8 reflective wavelength and normalized difference vegetation index continuity. Remote Sens. Environ. 2016, 185, 57–70.

34. Goodwin, N.R.; Magnussen, S.; Coops, N.C.; Wulder, M.A. Curve fitting of time-series Landsat imagery for characterizing a mountain pine beetle infestation. Int. J. Remote Sens. 2010, 31, 3263–3271.

35. Gandullo, J.M.; Sánchez Palomares, O. Estaciones Ecológicas de Los Pinares Españoles; ICONA: Madrid, Spain, 1994; ISBN 84-8014-116-6.

36. Aguirre-Salado, C.A.; Treviño-Garza, E.J.; Aguirre-Calderón, O.A.; Jiménez-Pérez, J.; González-Tagle, M.A.; Valdéz-Lazalde, J.R.; Sánchez-Díaz, G.; Haapanen, R.; Aguirre-Salado, A.I.; Miranda-Aragón, L. Mapping aboveground biomass by integrating geospatial and forest inventory data through a k-nearest neighbor strategy in North Central Mexico. J. Arid Land 2014, 6, 80–96.

37. R Core Team. R: A Language and Environment for Statistical Computing; R Foundation for Statistical Computing: Vienna, Austria, 2018.

38. Naimi, B.; Hamm, N.A.S.; Groen, T.A.; Skidmore, A.K.; Toxopeus, A.G. Where is positional uncertainty a problem for species distribution modelling? Ecography 2014, 37, 191–203.

39. Crookston, N.L.; Finley, A.O. yaImpute: An R package for kNN imputation. J. Stat. Softw. 2008, 23, 1–16.

40. Packalén, P.; Maltamo, M. The k-MSN method for the prediction of species-specific stand attributes using airborne laser scanning and aerial photographs. Remote Sens. Environ. 2007, 109, 328–341.

41. Hijmans, R.J.; van Etten, J.; Cheng, J.; Mattiuzzi, M.; Sumner, M.; Greenberg, J.A.; Lamigueiro, O.P.; Bevan, A.; Racine, E.B.; Shortridge, A. Package ‘raster. R Package 2019.

42. Dukes, J.S.; Pontius, J.; Orwig, D.; Garnas, J.R.; Rodgers, V.L.; Brazee, N.; Cooke, B.; Theoharides, K.A.; Stange, E.E.; Harrington, R.; et al. Responses of insect pests, pathogens, and invasive plant species to climate change in the forests of northeastern North America: What can we predict? This article is one of a selection of papers from NE Forests 2100: A Synthesis of Climate Change Impacts on Forests of the Northeastern US and Eastern Canada. Can. J. For. Res. 2009, 39, 231–248.

43. Camarero, J.J.; Gazol, A.; Sangüesa-Barreda, G.; Oliva, J.; Vicente-Serrano, S.M. To die or not to die: Early warnings of tree dieback in response to a severe drought. J. Ecol. 2015, 103, 44–57.

44. Marcone, C.; Franco-Lara, L.; Toševski, I. Major phytoplasma diseases of forest and urban trees. Phytoplasmas Plant Pathog. Bact.-I 2018, 287–312.

45. Ježić, M.; Poljak, I.; Šafarić, B.; Idžojtić, M.; Ćurković-Perica, M. ‘Candidatus Phytoplasma pini’ in pine species in Croatia. J. Plant Dis. Prot. 2013, 120, 160–163.

46. Senf, C.; Seidl, R.; Hostert, P. Remote sensing of forest insect disturbances: Current state and future directions. Int. J. Appl. Earth Obs. Geoinf. ITC J. 2017, 60, 49–60.

47. Klein, T.; Cohen, S.; Yakir, D. Hydraulic adjustments underlying drought resistance of Pinus halepensis. Tree Physiol. 2011, 31, 637–648.

48. Vicente-Serrano, S.M.; Lasanta, T.; Gracia, C. Aridification determines changes in forest growth in Pinus halepensis forests under semiarid Mediterranean climate conditions. Agric. For. Meteorol. 2010, 150, 614–628.

49. Bjorkman, C.; Niemela, P. Climate Change and Insect Pests; CABI: Oxfordshire, UK, 2015; Volume 8, ISBN 1-78064-378-0.

50. Weintraub, P.G.; Beanland, L. Insect vectors of phytoplasmas. Annu. Rev. Entomol. 2006, 51, 91–111.