Global Assessment of Climate-Driven Susceptibility to South American Leaf Blight of Rubber Using Emerging Hot Spot Analysis and Gridded Historical Daily Data

South American leaf blight (SALB) of Para rubber trees (Hevea brasiliensis Muell. Arg.) is a serious fungal disease that hinders rubber production in the Americas and raises concerns over the future of rubber cultivation in Asia and Africa. The existing evidence of the influence of weather conditions on SALB outbreaks in Brazil has motivated a number of assessment studies seeking to produce risk maps that illustrate this relationship. Subjects with dynamic and cyclical spatiotemporal features need to embody sufficiently fine spatial resolution and temporal granulation for both input data and outputs in order to be able to reveal the desired patterns. Here, we apply emerging hot spot analysis to three decades of gridded daily precipitation and surface relative humidity data to depict their temporal and geographical patterns in relation to the occurrence of weather conditions that may lead to the emergence of SALB. Inferential improvements through improved handling of the uncertainties and fine-scaled temporal breakdown of the analysis have been achieved in this study. We have overlaid maps of the potential distribution of rubber plantations with the resulting dynamic and static maps of the SALB hot spot analysis to highlight regions of distinctly high and low climatic susceptibility for the emergence of SALB. Our findings highlight the extent of low-risk areas that exist within the rubber growing areas outside of the 10◦ equatorial belt.


Introduction
Natural rubber is an essential raw material for the production of a wide range of goods and also the provision of numerous services (transportation above all).The Para rubber tree (Hevea brasiliensis Muell.Arg.) remains the unrivaled source for natural rubber.Although rubber trees originate from the Amazon basin, the Americas contribute only a fraction of the global rubber production (4.3% of 14.2 million tons in 2017) [1].This is due to South American leaf blight (SALB), the most serious disease affecting rubber trees.It is caused by the fungus Pseudocercospora ulei P. Henn.[2], which induces lesions on young leaves, leading to defoliation that, in multiple cycles, weakens and eventually destroys the trees.SALB has so far been confined to the Americas, but its threat haunts Asian and African rubber production.The established high-yielding clones that are the source for most of the rubber produced worldwide have a very narrow genetic base [3,4].They originate from a region in Para state of Brazil where the natural resistance to SALB among wild rubber trees is the lowest compared to populations growing in other areas within the endemic range of H. brasiliensis [5][6][7].The genetic erosion associated with decades of selection and breeding in a SALB-free environment has also contributed to the vulnerability of the Asian clones [8][9][10].By contrast, the evolutionary potential of P. ulei has been observed to be vigorous [11].The ongoing genetic improvement programs aiming to develop genotypes with stable resistance to SALB have yet to overcome this challenge [12].The destructive capacity of SALB, paired with the importance of its target crop [13] and the lack of effective and both environmentally and economically [14,15] acceptable countermeasures have earned P. ulei the status of a latent biological agent [16,17].Asian rubber-producing countries have implemented various measures to prevent the introduction of SALB and have developed contingency plans for its control in the case of an outbreak [15,[18][19][20].Regional conferences and workshops backed by the FAO (Asia and Pacific Plant Protection Commission) have aimed to increase knowledge sharing and action harmonization among the participating countries.
In the 1970s and 1980s, a series of programs encouraging natural rubber production (Programa de Incentivo à Produção de Borracha Vegetal-PROBOR) was conducted by the Brazilian government to increase national self-sufficiency, and although not successful, revealed that rubber trees could be grown in some areas of Brazil (mainly in some parts of São Paulo State) without being affected by SALB [21].These "escape zones" are colder and experience an annual dry period that coincides with the refoliation period of the rubber trees [22].
The 'disease triangle + time' conceptual model in plant pathology [23,24] describes disease emergence as the result of the convergence of a susceptible host and a pathogen under favorable environmental conditions over a sufficient period of time.The absence of any of those elements would restrain disease development.Humidity, temperature, and light conditions have been reported to play a role in creating the conditions for SALB infection by influencing release, viability, and germination of spores and facilitating penetration of host tissues [15,[25][26][27][28][29][30]].Among them, long lasting humid conditions leading to several hours of continuous "leaf wetness" is the point on which consensus among old and new literature exists.
Mature rubber trees go through an annual defoliation and refoliation cycle as opposed to seedlings, which within the first four to five years produce multiple leaf flushes throughout the growing season.Young rubber leaves are susceptible to SALB until up to three weeks after bud-break [31][32][33] when they reach their physiological and structural maturity [34,35].
A few studies based in geographic information systems (GIS) exist that focus on mapping SALB "escape zones", but on a relatively small scale [22,[36][37][38].Roy et al. [39] expanded the covered geographical evaluation span to a global scale and took the effects of different phenological behavior of mature and immature trees into consideration.However, the limitation in this study is the time steps, which are monthly, and considering the pace at which the refoliation process and SALB infection evolve, lacks the requisite detail.The annual and monthly variables and their corresponding thresholds used in comparable studies are practically used as proxies attempting to capture the effects of the underlying daily variations.However, thanks to improved access to spatiotemporal historical data with finer temporal granulation (see Section 2.1) and an increase in computational power, the use of such proxies is more and more avoidable.Extraction of trends from weather records should account for and present the inherent uncertainties for highly dynamic variables, such as precipitation.
In this study, we applied emerging hot spot analysis, which is based on the Getis-Ord Gi* statistic [40,41], to three decades of gridded daily precipitation and surface relative humidity data in order to reveal temporal (in daily steps over the yearly cycle) and geographical patterns that link weather conditions to the likelihood of SALB emergence.

Data
We used the CPC (Climate Prediction Center) global unified gauge-based analysis of daily precipitation data [42] and the NCEP (National Center for Environmental Prediction) global reanalysis daily surface relative humidity data [43], which are both available at the website of the Physical Sciences Division of the Earth System Research Laboratory under the Office of Oceanic and Atmospheric Research-National Oceanic and Atmospheric Administration (NOAA/OAR/ESRL PSD).These two data sources have been evaluated and used for long-term trend analysis of humidity and precipitation with reasonable results in different parts of the world [44][45][46][47][48].The CPC precipitation data and the NCEP relative humidity data have spatial resolutions of 0.5 and 2.5 decimal degrees, respectively.We chose a 30-year period from 1988 to 2017 as the time frame for our input data.Based on the global climatic suitability map for rubber cultivation (Figure 1), we limited the geographical coverage of this study to the 30 • north and south latitudes in which 0.5 decimal degrees range from about 55.6 km at the equator to 47.3 km at 30 • .Global distribution of climatic suitability for rubber cultivation.This map is based on four temperature and precipitation criteria: annual mean temperature, intra-annual temperature distribution, annual precipitation, and intra-annual precipitation distribution [22] applied to Worldclim V2 [49] gridded climatic data.The process for generation of this map is described in detail in Golbon et al. [50].

Methods
Our aim here is to explore the spatiotemporal patterns for high relative humidity and rainfall, which, when combined, represent the long-lasting hours of leaf wetness, facilitating SALB establishment.Daily precipitation and relative humidity records in the continuous form were extracted from the NetCDF (Network Common Data Form) files.Surface relative humidity layers were resampled to match the spatial resolution of the precipitation data.Using conditional expressions applying thresholds of 1 mm for precipitation (representing precipitation incidence regardless of the quantity) and 90% for relative humidity (standing for close to saturation air humidity [39]), we converted the continuous data to binary layers (0 for values lower than the threshold and 1 for otherwise).These binary raster files were converted to point shapefiles.The points corresponding to waterbodies and areas with latitudes beyond 30 • N and 30 • S were eliminated.We added a new date-type field to the attribute table of each point shapefile and populated it with the corresponding date value.Point shapefiles affiliated with each day of the year (e.g., all the 30 files related to Jan 01 from 1988 to 2017) were merged, and the resulting shapefile was projected using the Hammer-Aitoff projection, which is metric and global.The Emerging Hot Spot Analysis (EHSA) tool of ArcGIS does not process the usual raster and vector inputs directly.Point shapefiles containing the temporal information (described above) need to be aggregated into NetCDF data structures called space-time cubes, which in addition to the coordinates, stores the time dimension information.This three-dimensional data structure has spatial units (bins) that can be in the form of either squares or hexagons.The projected point shapefiles were used to create space-time cubes containing the temporal and spatial information on events of precipitation and high surface relative humidity.We chose the hexagonal form for the space-time bins and, based on the resolution of the precipitation data (see Section 2.1), fixed the neighborhood distance (the height of the hexagons) to 90 km.The space-time cubes were used as input for the EHSA.
EHSA uses the Getis-Ord Gi* statistic [40,41] to evaluate the spatiotemporal patterns of occurrence/absence of an event of interest (precipitation and high relative humidity in the current case), verifying them versus the probability of the observations being outcomes of random processes [40,51,52].EHSA returns three possible main categories for each space-time bin: cold spot, hot spot, or no pattern.In the context of this study, area units (hexagons) defined as hot spots are rendered to be associated with conditions favorable to the pathogen activity (incidence of rainfall or high relative humidity) as opposed to cold spots, for which the investigated criterion was detected not to be within the favorable range for the pathogen.The 'no pattern' category allows for the indeterminate conditions in which the data at hand supports neither a hot nor a cold spot category.The cold and hot spots were each divided into eight subcategories differing in their observed temporal stability and if applicable, the shifting direction detected for them over time.However, for the purposes of this study, we focused on the hot and cold spots with significant trends (Mann-Kendall trend test z-scores with p-values smaller than 0.05) regardless of their subcategories in the subsequent steps and referred to them with acronyms "HS" (hot spot) and "CS" (cold spot), respectively.We handled the areas for which no pattern was detected and those with statistically nonsignificant trends equally and labeled them with "NS" (no pattern/nonsignificant trend).The EHSA output attributes were simplified from the original 17 possible outcomes to these three simplified categories.In order to homogenize the day by day changes in the geographical delineation of the HS/CS/NS hexagons, we extracted the majority outcome for each hexagon over a moving seven-day time window centered on the target day of the year (e.g., 28 Dec to 04 Jan for 01 Jan).
The day-to-day maps for both relative humidity and precipitation were date-tagged, laid out and exported as image files.We used Blender 2.79 to render an animated map with these 365 frames to be played as a 15 second long video.We further simplified the dynamic map in two steps.First, by adding a new field to the attribute table of each shapefile and merging the HS/CS/NS information of both criteria from which six possible classes ensued (HS + HS, HS + NS, NS + NS, CS + NS, CS + CS, and HS + CS).Finally, we added another field to the attribute tables and populated them with two possible classes: hexagons for which at least one of the two criteria was a significant cold spot (=1) and otherwise (=0) as the alternative.Finally, the resulting binary values for each hexagon were summed up over the whole year and for the refoliation periods to produce two static summary maps.In order to secure the coverage of the refoliation periods, we applied the time window of Feb-Apr for the Northern Hemisphere and Aug-Oct for the Southern Hemisphere.We have also created keyhole markup language zipped (KMZ) files corresponding to all of the daily maps from which the three mentioned videos were produced.

Results
The key outputs of this study are the daily spatial classifications presented as three animated maps (Videos S1 to S3), which require some direct engagement by the users (navigating through the time frames for an area of interest) to reveal the desired information.Video S1 presents the detected trends for both relative humidity and precipitation, while Video S2 simplifies that content by reflecting only the six possible class combinations regardless of their origin (i.e., which class belongs to which criterion).Video S3 focuses on the significant cold spots and combines the hot spots with the no pattern/not significant class.Although further compression of the outputs into single frame maps summarizing the patterns occurring over longer time periods may seem to be in conflict (by forcing potentially arbitrary time frames and risking overgeneralization) with the goals of such a study, it is prudent to describe the extremes observed over the complete year cycle and the approximate leaf flushing period (Figure 2).Areas for which no significant cold spots (for neither of the two criteria) were detected across the yearly cycle were confined to the Asian and American tropics between 9 • N and 8 • S (Figure 2a).The other extreme outcome-the entire yearly cycle covered by a significant cold spot for at least one of the two criteria-was mostly returned for Asian, Pacific, and Caribbean islands and peninsular relatively narrow land strips.
Limiting the time window of interest to the approximate refoliation period for mature trees (Figure 2b), we found an increase in the area for which no significant cold spots were detected in the American tropics from the equator to the 8 • N latitude.Within the same temporal frame, this pattern emerged around the African equatorial line from longitudes 8 • E to 20 • E. The outcome for the aforementioned insular and peninsular areas was very similar to findings on the annual cycle.In addition to what was indicated for the year-round summary, the following areas were found to be characterized throughout the refoliation period by at least one of the two criteria as a significant cold spot: in Asia, coastlines of the two Indian states of Kerala and Karnataka (8 • N to 14 • N) and the eastern Indian coasts from 11 • N to 13 • N and 18 • N to 21 • N; in Africa, the area enclosed within 11 • N, 15 • W and 31 • E, 4 • N and at the eastern African coastlines from 5 • S to 8 • S (Tanzania) and 17-18 • S (Mozambique); and in South America, northern Venezuela, Brazilian coastline from 6 • S to 18 • S and 42 • W to 50 • W stretching deep into Maranhão and Pará states and a large cluster in the southern parts of Mato Grosso state.In general, the share of the two extremes (all days vs. no days associated with a significant cold spot) from the total area climatically conducive to rubber cultivation is larger for the refoliation period than the annual cycle, reflected by the higher contrasts in panel b than a of Figure 2 and the transition from potential risk zones to safe zones occurs in shorter distances.Days for which at least one of the two criteria considered in this study was detected to be a significant cold spot from the complete year cycle (a) and the approximate refoliation period (b) (Feb-Apr for the Northern Hemisphere and Aug-Oct in the Southern Hemisphere) are reflected in a proportional form.In this figure, similar to Video S3, we have pooled hot spots and the uncertain areas distinguishing them from cold spots which form an 'unsafe/uncertain' vs. 'safe' dichotomy.Therefore, 0% cold spot means 100% hot spot/uncertain.The non-hatched area reflects the zones climatically conducive to rubber cultivation (for more details see Figure 1).

Discussion
The climatic risk analysis of SALB by Roy et al. [39] is the only publication comparable to this study in regards to the time dimension during data processing and inference, the similarities in the covered geographical span, and the similarly rough spatial resolution of the input data.That study diverges from this work in three major points: exclusion of temperature variables from analysis, data processing method (EHSA vs. rule-based classification), and the level of temporal granulation (daily vs. monthly time slices).The temperature ranges reported in scientific literature for P. ulei [25,53,54] focus mainly on the optimal conditions for the pathogen, which cannot be used as impeding thresholds.Moreover, low temperature appears not to be a limiting factor for P. ulei within the geographical range that is climatically suitable/bearable for H. brasiliensis [26].Figure S1 in the Supplementary Materials demonstrates how the spatial limits set by relative humidity and precipitation are completely nested within the space outlined by the temperature limits.Therefore, temperature turns to a superfluous criterion.There are precedents for the temporary loss of immunity to SALB in some established Brazilian 'escape zones' due to unusually humid weather conditions [21], but cases relating to temperature anomalies have not been reported.
We chose to use EHSA, which not only accommodates temporal trends but also can allow for uncertainties and avoid potentially overgeneralized inferences based on vague dichotomies.Technical barriers, such as insufficient access to data and computational power, which are often the reasons for resorting to using proxy variables, are being progressively alleviated thanks to recent advances.
As opposed to small-scale and local studies, investigations with global coverage help provide a bigger picture.This broad overview comes at the cost of the potential loss of details compared to the reality, such as missed smaller but not necessarily insignificant pockets of land with conditions contrasting the dominant neighboring land use and shifts in some of the detected borderlines.The broad overview mentioned above is helpful in evaluating the plausibility of the findings and the performance consistency of the data and methods in a general glimpse and also as an informal cross-validation.For instance, persistent cold spots assigned to the Caribbean and small Asian islands (see Figure 1 and Videos S1 to S3) are dubious considering the stagnation of the class attributed to them over the year cycle and in relation to the conditions reflected in their immediate vicinity.This situation is probably caused by spatial isolation (shortage of neighboring points) of the corresponding space-time bins.
Ideally, a map that encapsulates the trends for the approximate annual refoliation period should restrict the summarized values assigned to each grid cell to the dates specifically relevant for the corresponding location.The surprisingly scarce and to some extent speculative literature that exists in regards to the processes and triggers of the phenological behavior of Para rubber [55][56][57] does not aggregate to a solid base to link refoliation to the geography through the natural annual cycles.As an alternative, metadata gathered from studies that report both refoliation dates and the corresponding geographical coordinates from different locations should serve as an alternative for producing a time map for refoliation.The existing literature which report both refoliation dates and coordinates [58][59][60][61][62][63][64][65][66][67] are not abundant enough and the indicated refoliation dates are often formulated in too much a loose range to be applicable to our purpose.As a result, the map which summarizes the situation over the refoliation period (Figure 2b) inevitably suffers from overgeneralization of the time frame.Nevertheless, this problem is tackled using the dynamic maps (Videos S1 to S3).This requires users to navigate through the dates (video frames) of interest for a given area.A media player that permits the frame by frame viewing of the video should be used for that purpose.We suggest GOM media player (available at www.gomlab.com),which permits frame by frame navigation by pressing the 'f' key on the keyboard.
Roy et al. [39] reported all mature trees in Indian rubber growing areas to be safe from SALB.Our findings support this assessment for the Western Ghats, where cold spots disappear in mid-May and reappear in November, and for the Odisha State where the low humidity period is even longer, stretching from October to June.However, the identified trends in the northeastern region of India are associated with considerable uncertainty.Moving further to the east, Roy et al. [39] describe most of the rubber growing areas in China, Thailand, Cambodia, Laos, and Vietnam as low-risk, and parts of Philippines, Malaysia, and Indonesia as high-risk zones.We reached comparable outcomes for these areas.For the continental parts of Southeast Asia, however, the time gap between the end of the refoliation period and the start of the humid season tends to be narrow.Roy et al. [39] applied a blanket refoliation time window of February-March for their assessment of SALB risk for mature rubber trees in Africa.In areas south of the equator, this is not consistent with their refoliation rule of thumb, March-May for the Northern Hemisphere and September-November for the Southern Hemisphere.With using the dual time window, our findings suggest that the conditions during the refoliation season over large areas north of 5 • N from Guinea to the Central African Republic and east of 30 • E from Mozambique to Kenya (Figure 2b) tend to not support leaf wetness.
The thresholds applied in this study represent conditions leading to long hours of leaf wetness through the coincidence of rainfall and high relative humidity.However, other conditions causing leaf wetness (e.g.dew in place of rainfall) are also possible.Despite the fact that by the three levels of simplification applied to the dynamic maps some flexibility for more stringent examination of the outputs has been preserved, the thresholds applied to the raw data are fixed.
For studies relating to subjects concerning alien species, such as this one, the following seven measures suggested by Groom et al. [68] should be followed: (1) creation of data management plans, (2) increasing interoperability of information sources, (3) documentation of data through metadata, (4) formatting data using existing standards, (5) adoption of controlled vocabularies, (6) increasing data availability, and ( 7) ensuring long-term data preservation.We have achieved points ( 2) and ( 4) by providing the KMZ outputs of our study to the public, and points (3), ( 5), (6), and ( 7) by archiving these files on zenodo.org(https://doi.org/10.5281/zenodo.2576857)data repository.These KMZ files can be used by any interested party regardless of their expertise, examined by individuals of no familiarity with geographic information systems (dragging and dropping in GoogleEarth™), or exploited by GIS-savvy researchers for further developments.

Conclusions
Incorporation of daily time steps in the risk evaluation is a tedious process but achievable.The uncertainties rising from the natural overtime variabilities in observations can be accounted for using EHSA.For rubber trees past their multiple-leaf-flush life stage, several areas safe from SALB (regarding the leaf wetness during the refoliation period) appear to exist in the areas outside of the 10 • equatorial belt.

Supplementary Materials:
The following are available online at http://www.mdpi.com/1999-4907/10/3/203/s1, Figure S1: The de facto redundancy of monthly mean temperature in detection of areas climatically suitable for the formation of South American leaf blight (SALB); Video S1: Trends for humidity and precipitation in relation to outbreak risk of South American leaf blight; Video S2: Simplified trends for humidity and precipitation in relation to outbreak risk of South American leaf blight; Video S3: Significant risk cold spots for South American leaf blight outbreak overlaid with climatically optimal and suboptimal areas for rubber cultivation.
Figure1.Global distribution of climatic suitability for rubber cultivation.This map is based on four temperature and precipitation criteria: annual mean temperature, intra-annual temperature distribution, annual precipitation, and intra-annual precipitation distribution[22] applied to Worldclim V2[49] gridded climatic data.The process for generation of this map is described in detail in Golbon et al.[50].

50 Figure 2 .
Figure 2.Summary maps of the cold spot occurrence.Days for which at least one of the two criteria considered in this study was detected to be a significant cold spot from the complete year cycle (a) and the approximate refoliation period (b) (Feb-Apr for the Northern Hemisphere and Aug-Oct in the Southern Hemisphere) are reflected in a proportional form.In this figure, similar to Video S3, we have pooled hot spots and the uncertain areas distinguishing them from cold spots which form an 'unsafe/uncertain' vs. 'safe' dichotomy.Therefore, 0% cold spot means 100% hot spot/uncertain.The non-hatched area reflects the zones climatically conducive to rubber cultivation (for more details see Figure1).