The red palm weevil (RPW), Rhynchophorus ferrugineus (Olivier) (Coleoptera: Curculionidae), is an economically important pest of various palm species that has invaded all continents (Kaakeh 2005; Faleiro 2006; Ju et al. 2010; NAPPO 2010). It appears to be native to South and Southeast Asia. The pest was first described as a deadly coconut pest in northern India (Lefroy 1906) then later reported on date palms (Madan 1917). By the mid 1980s, the pest had spread to the Middle East (Gomez & Ferry 1999; Abraham et al. 2000; Faleiro 2006), and then it moved more rapidly into northern Africa by 1992 and southern Europe by 1994, eventually reaching North America in 2009 (Cox 1993; Abozuhairah et al. 1996; Barranco et al. 1996; Faghih 1996; El-Ezaby 1997; Kehat 1999; Murphy & Briscoe 1999; Quin et al. 2002; EPPO 2007, 2008, 2009, 2010; Kaakeh 2005; Kontodimas et al. 2006; Al-Eryan 2009; Alhudaib 2009; Borchert 2009; Bertone et al. 2010; Ju et al. 2010; Roda et al. 2011). The pest is reported from 30 palm species (Murphy & Briscoe 1999; Faleiro 2006; Kontodimas et al. 2006; EPPO 2008; Bertone et al. 2010; Dembilio & Jacas, 2012) used in agriculture and landscaping. The larva feeds within the trunk and by the time the first symptoms are detected, the attacked palm is already seriously injured and the population of the pest reaches already high levels (El-Ezaby 1997; Faleiro 2006). This cryptic phase protects the pest from direct harsh external climatic conditions and therefore enables it to breed in a wide range of climates (Murphy & Briscoe 1999), and at the same time enables easy spread through vegetative planting materials. Huge shipments of planting material from one country to another one have therefore contributed tremendously to the rapid spread of the pest (Faleiro 2006).

A single RPW female can lay 58–760 eggs in its lifetime (Avand-Faghih 1996; Abraham et al. 2002; Kaakeh 2005; Faleiro 2006, Prabhu & Patil 2009). The adult is a strong flier that can move >900 m in a single flight, and as much as 7 km in 3–5 d (Abbas et al. 2006). The larvae (grubs), which are the destructive stage of the pest, live 25–105 d before pupating. Various studies of effects of temperature on its development (Salama et al. 2002; Martin & Cabello 2006; Dembilio & Jacas 2011; Li et al. 2010) have yielded diverse results. The studies found minimum temperature tolerances ranging from as high as 17.4 °C to as low as -2.3 °C for RPW (Salama et al. 2002; Martin & Cabello 2006; Dembilio et al. 2010 and Li et al. 2010).

The cryptic life of the pest makes its control difficult, and eradication of the pest has not been achieved in any of the invaded areas. The current management approach involves an integrated pest management (IPM) composed of monitoring, mass trapping, insecticide application, and early detection (Murphy & Briscoe 1999; Faleiro 2006). Few natural enemies have been found associated in the field and although some were found effective under laboratory conditions, none has so far led to successful control of the pest per se (Murphy & Briscoe 1999; Faleiro 2006; Shahina et al. 2009; Dembilio et al. 2010; Dembilio & Jacas 2011).

Invasive species are an issue of great concern globally, particularly in light of the ever-increasing scale of human movement and trade globalization (Levine & D'Antonio 2003; Hulme 2009). Given the economic impact and continuing spread of RPW, there is a considerable need to being able to anticipate new areas of the pest's invasion. Proactively identifying locations suitable for its establishment may enable decision-makers and agricultural and environmental protection officers to initiate preventative measures or rapid responses in timely manner. A recently developed approach to predicting species distribution is ecological niche modeling (“ENM”; Peterson 2003) where a suite of techniques are used to estimate the species' environmental requirements in broad, coarse-resolution dimensions (Soberén 2007, 2010). Once estimated and the estimate evaluated for predictive ability, the niche model can then be projected onto other regions to identify areas matching inhabited areas that may represent potentially new areas of invasion (Peterson 2003; Peterson & Vieglais 2001; Sutherst & Maywald 2005; Fiaboe et al. 2006).

Here, we use ENM approaches to explore the likely geographic potential of this invasive species. Because the performance can differ among methods (Elith et al. 2006;Peterson et al. 2008), we explored the predictive ability of 2 commonly used algorithms: the desktop version of the Genetic Algorithm for Rule-set Prediction (GARP; Stockwell & Noble 1992) and MaxEnt (Phillips et al. 2006). We subdivided the present known range of the species into calibration (training) and evaluation (testing) areas to test predictions rigorously across unsampled landscapes, to assure that the models have predictive power. These tests establish that the species follows environmental rules that are consistent across multiple continents that are likely very different in species composition and environments (Jiménez-Valverde et al. 2011). With validation of the predictive ability of the model, we projected the ENMs worldwide to identify areas of potential distribution and possible invasion for the species.

METHODS

Known Distribution of RPW

The initial occurrence data set included 132 localities at which red palm weevils are known to have occurred around the world. These points were compiled from various publications, personal communications, and our own research (Appendix 1; Fig. 1). Because geographic coordinates of occurrence points were reported in the literature for only a few locations, we generated coordinates via reference to the directory of world named places in the Global Gazetteer version 2.2 ( http://www.fallingrain.com/world; see summary of known occurrences in Fig. 1). The single, recent occurrence in the United States was eliminated from the data set to provide a level of independence, as the distributional potential of the species in the United States was the immediate impetus for this study.

Fig. 1.

Example (replicate 1) of spatially stratified predictions into 2 regions in the second and fourth quintiles of longitude in the known present distribution of red palm weevils (summary of overall sampling regime shown in top panel). Occurrences are shown as Xs. Model predictions are shown as ramps from white (= unsuitable) to dark orange (highly suitable). A = World View of Occurrences, 200 km Buffer, and Testing Areas; B = Western Testing Region — Iteration 1; and C = Eastern Testing Region, Iteration 2.

Originally, these occurrences were not distributed uniformly across world landscapes, but rather were highly clumped in their distribution (e.g., large clusters of points in southern India). Hence, to avoid pseudoreplication of local environments owing to artificial clustering of occurrence sites, we reduced the raw set of localities to 38 spatial clusters, each separated from each other by >100 km. We then represented each of these clusters once in each of 5 replicate data sets, choosing random sets of representatives of each cluster in each replicate model.

Niche Modeling and Calibration

Ecological niche models are ideally fitted within the area that has been accessible to the species over time periods relevant to its distributional history (Barve et al. 2011).Contrasts are made between environments associated with known occurrences and those associated with sites at which the species is not known to occur. We based the models on climatic features that are related to species' natural history, particularly parameters related to heat, cold and water stress. From the broader suite of “bioclimatic” parameters available worldwide (Hijmans et al. 2005), we chose 7 that are relatively uncorrelated globally (Jimenez-Valverde et al. 2009): annual mean temperature, mean diurnal temperature range, maximum temperature of warmest quarter, minimum temperature of the coldest quarter, annual precipitation, precipitation of wettest quarter and precipitation of driest quarter (Beaumont et al. 2005). In light of the spatial precision of the distributional data available, where we had few or no data that were finer than ∼5km in terms of their known spatial distribution, we chose 2.5' spatial resolution as most appropriate for our analyses to avoid over interpreting the data.

A framework for understanding distributions of species is termed the “BAM” diagram (Soberón & Peterson 2005), in which the species biotic, abiotic, and mobility constraints are estimated to the extent possible; the distribution of the species is in essence - the intersection of the 3 areas. Of particular relevance is the mobility constraint (the area termed “M”), which is the area that the species has sampled through time, and within which the species can be assumed to have colonized all sites presenting suitable conditions (Barve et al. 2011). The invasive nature of RPW makes a dispersal-focused M definition appropriate (Barve et al. 2011). As a consequence, we buffered all occurrence points by 200 km to create an arena for modeling, as an approximation of the area that has been accessible to the species over the recent past. To reduce the degree to which spatial auto-correlation might compromise model testing, we used a spatial subsetting exercise. Specifically, we divided the 38 occurrence areas into quintiles by longitude, each of which held 7–8 occurrence areas. We used the first, third, and fifth quintiles to train models and the other 2 areas to test predictions (see Fig. 2); the arena for model evaluation was thus the union of these 2 areas, but only within the hypothesis of M described above.

We estimated ecological niches using the 2 niche modeling algorithms that are perhaps the 2 that have seen the most use in the literature, MaxEnt (Phillips et al. 2006) and GARP (Stockwell & Peters 1999). MaxEnt is a method developed to estimate ecological niches of species based only on presence data, although the broader ‘background’of conditions across the study area is used in the analysis. The information available generally takes the form of a set of real-number-valued environmental variables, called “features,” and distributions are fitted under the constraint that expected values of each feature should match the empirical average (average value for a set of sample points taken from the target distribution). MaxEnt thus attempts to estimate the probability distribution for the occurrence of species as the “maximum entropy” distribution. The result is an approximation to a uniform probability distribution, subject to the constraints imposed by the environmental conditions associated with known occurrences of the species in question (Phillips et al. 2006).

MaxEnt is relatively robust to small sample sizes, but sites sampled must represent the environmental diversity of the species and the study area for models to be robust (Pearson et al. 2007; Wisz et al. 2008). A real-number suitability value is assigned to each pixel, which can vary from 0 (no suitability) to 1 (complete suitability). To avoid overfitting (i.e., avoiding predictions that fit well to training data but have little generality), raw continuous predictions were converted to binary formats by means of a thresholding step explained below. Maxent version 3.3.1 with the random seed option, logistic output options, and sample-size-dependent feature choice were used. However, 50% of occurrence data were reserved for testing, and a random seed was used to assure distinct runs in subsequent tests. MaxEnt output was imported into ArcGIS 10 as floating point grids. In Arc 10, we multiplied these raw grids by 1000, and truncated them to create integer grids.

GARP searches complex solution spaces using a genetic algorithm. Within GARP processing, input occurrence data are divided randomly into 3 subsets: training data (25%; for model rule development), intrinsic testing data (25%; for intrinsic evaluation and tuning of model rules) and extrinsic testing data (50%; for evaluation of model quality and filtering among replicate models, see below). Spatial predictions of presence versus absence can include 2 types of error: false negatives (areas of actual presence predicted absent) and false positives (areas of actual absence predicted present; Fielding & Bell 1997); rule performance in terms of overall error is evaluated via the intrinsic testing data set. Changes in predictive accuracy from one iteration to the next are used to evaluate whether particular rules should be incorporated into the model or not, and the algorithm runs either 1000 iterations or until convergence (Stockwell & Peters 1999). The final rule set is then used to query the environmental data sets across the study region to identify areas fitting the rule-set prediction, producing a hypothesis of the potential geographic distribution of RPW (Soberón & Peterson 2005).

Fig. 2.

Summary of global projections of ecological niche models trained based on the known occurrences of red palm weevils globally (except for the California occurrence, which was omitted from analyses). A = predicted suitability, on a ramp from white (unsuitable) to red (highly suitable). B = the degree of novelty of the environments represented, with blue indicating environments closely similar to the points of known occurrence, and red indicating environments that are widely different; model projections into regions at the latter end of this novelty scale are suspect.

Since GARP processing includes several random-walk components, each replicate model run produces distinct results, representing alternative solutions to the optimization challenge. Consequently, following recommended consensus approaches (Anderson et al. 2003), we developed 100 replicate versions of each model. We filtered these replicates based on their error characteristics to emphasize the overriding importance of omission error (as opposed to commission error), retaining the 20 models showing the lowest false-negative rates (= percentage of independent testing points falling in areas not predicted to be suitable), and then retaining the 10 (of the 20) closest to the median of the proportional area predicted present among models, an index of false-positive error rates (Anderson et al. 2003). A consensus of these ‘best subset’ models was then developed by summing values for each pixel in the map to produce final predictions of potential distributions with 11 thresholds (i.e., integers from 0 to 10).

Model evaluation.—Once predictions were developed with MaxEnt and GARP, we reduced the predictions to leave only the testing region using the Extract by Mask feature of ArcGIS, version 10. We exported the attributes table associated with this raster as a summary of proportional areas predicted as suitable at each predictive level. We extracted raster values to the testing points to assign predictive levels to each point, and by extension establish omission error rates to each predictive level of the model. We used a partial ROC approach that allows reweighting of error components in a ROC framework, emphasizing the dominant role of omission over commission error in evaluations of model quality (Peterson et al. 2008). This method is designed around a parameter E that estimates how much environmentaly significant positional error is likely present in the occurrence data—essentially how much omission error would be expected if the model estimated the habitable areas of the species perfected. In light of the nature of the data that we used in this study (“found” data), we used E = 0.1, a relatively high error rate, to define the portion of ROC space within which to evaluate model predictions. For each model test, we used direct count of ROC area under the curve (AUC) scores out of a bootstrapped resampling of 50% of available testing data as an estimate of the probability associated with the particular model prediction.

RESULTS

Our literature search identified 132 sites at which RPW is known to occur worldwide (Appendix 1), shown in Fig. 1. These points form the basis of all of our model development. Detailed evaluations of model predictions into independent testing regions indicated that MaxEnt models yielded predictions that were statistically significantly elevated above random expectations for all 5 random replicate subsamples (all P < 0.002; Fig. 1). GARP models failed to achieve statistical significance in 3 of 5 cases, in which predictions achieved probabilities of 0.07 < P < 0.10, outside of the range of statistical significance. As a consequence, we used only MaxEnt predictions in the remainder of the analyses in this study, and Max-Ent predictions were amply confirmed as having robust predictive power regarding the potential geographic distribution of this species, even in broad regions from which no occurrence data were available.

Global projections, effectively hypotheses of environmental suitability of landscapes based on environmental characteristics of known sites of occurrences, indicated a pantropical potential distribution for the species, ranging from East and Southeast Asia westward across the Indian Subcontinent to West and Central Africa and northern South America (Fig.2A). MaxEnt also identified areas at high northern latitudes as suitable, albeit only at moderate levels. However, the MESS maps (Fig.2B) indicate that these areas were remote in environmental space from the set of conditions under which models were calibrated, indicating that these extrapolations should be accorded little weight.

Touring around the world for potential RPW distributional areas, and bearing in mind that the models were calibrated using data from these same regions, in South and Southeastern Asia, model predictions covered the known distribution of the pest in each country where the pest has been reported. The model did, however, identify suitable areas in Nepal and Bhutan where the pest has not as yet been reported. In the Middle East, all of the known distribution was predicted, except for areas of Georgia and Iraq that are known to hold infestations. Areas of the Caucasus not currently known to hold the pest (Azerbaijan, Armenia), were predicted as suitable for establishment of populations of the species (Fig. 2).

In Europe, the model replicated all known distributional areas, but also extended farther north, to include portions of Belgium, Denmark, Estonia, Finland, Ireland, Latvia, Lithuania, Netherland, Norway, Poland, Sweden, and United Kingdom; the biological reality of these predictions is uncertain (Figs. 2 and 3). Across Africa, all known distributional areas were included in model predictions; the Sahara Desert was predicted as unsuitable, and the Sahel and southern Namibia and Botswana were identified as areas of relatively low suitability; however, all other Subsaharan African countries were predicted as holding conditions highly susceptible for this pest.

In the Caribbean islands, Aruba and Curaçao, where the pest is already established, were predicted as suitable; in addition, all Caribbean islands were predicted as suitable, suggesting considerable potential for spread in this region (Figs. 2 and 3). In South America, highly suitable zones for establishment were identified extending from Venezuela and Colombia south to Bolivia and northern Argentina. All Central American countries were found to present suitable conditions for red palm weevil establishment.

In North America, suitable areas were identified across the southern United States and Mexico. As mentioned above, and following the novel environment map (Fig. 2), the apparently suitable areas in northern Canada and Greenland are highly extrapolative and should not be considered as suitable for establishment. In the United States, the single occurrence point (Laguna Beach, California) was successfully predicted by the model (recall that this point was omitted from calibration datasets) (Fig. 3). Additional areas that were reconstructed as highly suitable for pest establishment included 16 counties in Florida, 4 in Louisiana, and coastal portions of 13 counties in California; overall, however, the potential for RPW establishment in the USA is limited (Fig. 3).

DISCUSSION

The niche modeling methods utilized herein provide a useful, independent view on the laboratory-based results regarding temperature contrasts reported by various authors (Salama et al. 2002; Martin & Cabello 2006; Dembilio & Jacas 2011; Li et al. 2010). Li et al. (2010), studying Chinese populations, reported a lowest temperature for successful development of 17.4 °C, with an accumulated temperature of 1,590 degree-days (DD) required. These authors also reported the pupal stage as the most resistant to cold, with a lower thermal threshold of 16.5 °C. Martín and Cabello (2006) reported similar trends in laboratory studies of Spanish populations, with thermal thresholds of 13 °C and 15 °C for the pupae and larvae, respectively, and 1,436 DD from larva to adult hatching. However, Dembilio & Jacas (2011), studying Spanish population in live palms in a greenhouse setting, reported mean monthly thermal thresholds as low as 4.5 °C for the second larval instar to pupa, with a much-lower total thermal constant of development of 989.4 DD from egg to adult. Dembilio et al. (2011) reported thermal thresholds of 15.45 and 13.95 °C for RPW oviposition and egg hatching. Salama et al. (2002) reported in Egypt a low thermal threshold of -2.3 ° C for the pupal stage of the pest. In Egypt, El Ezaby (1997) reported an upper temperature threshold of eggs for hatching at 40 °C.The niche models, provide a view that is quite independent, based on geographic and environmental range limits rather than on individual and population tolerances.

One complication to the niche model results, however, is that of extrapolating model predictions from the relatively restricted known distributional areas and associated estimate of M for the species to areas worldwide. The MESS maps implemented in MaxEnt (Elith et al. 2011) provide some insight into these areas of extrapolation, but also commit a logical error. The MESS maps summarize distance in environmental space to the known occurrence points; however, if our calibration area matches the M for the species (Barve et al. 2011), and if that area is the area that the species has “sampled” over its history (and colonized or not, given conditions manifested there), then the MESS maps should instead contrast global environments to those present across M. Otherwise, the result confuses the environmental limitations inherent in the ecological niche of a species with the environmental limitations of the input data and calibration process. Improved versions of MESS, and routines for their convenient estimation and implementation, are under development (J. Soberón, University of Kansas, personal communication).

Our models indicate large areas of the world that remain susceptible to RPW invasion, such as much of Subsaharan Africa, the Caribbean, and other areas. The model predicted successfully all knows distribution. In addition, the prediction for China is similar to the potential establishment obtained by Ju & Ajlan (2011) who used a phenology model approach for China. Although model specificity (i.e., avoidance of including “extra” areas in the prediction of suitable areas) is a concern, particularly in light of the possible overextension of model predictions in Europe, northern Canada and Greenland, we suspect that the invasion of the rest of the world by this species remains incomplete. The novel environment map is therefore needed for practical use of the predicted suitability model. In addition, the presence of host plant will also be a tremendous guide for use of the present suitability model in crop protection programs whether for survey and monitoring, quarantine or management of the pest. The natural host range of RPW covers primarily the palms (Arecaceae). A total of 32 plant species belonging to 3 families (Agavaceae, Arecaceae and Poaceae) were reported as suitable host (Murphy & Briscoe 1999; Faleiro 2006; Kontodimas et al. 2006; EPPO 2008; Malumphy & Molan 2009; Bertone et al. 2010; Dembilio & Jacas 2012). In the Arecaceae family, a total of 30 plant species were recorded and including: Areca catechu L., Arenga saccharifera Labill. ex DO, A. pinnata (Wurmb) Merr., Borassus flabellifer (Mart.) Warb., Borassus sp., Brahea armata S.Watson, Butia capitata (Mart.) Becc, Calamus merrillii, Caryota cumingii Lodd. Ex Mart, C. maxima, Cocos nucifera L., Corypha utan Lam. (= C. gebanga, C. elata), C. umbraculifer L., Elaeis guineensis Jacq., Livistona decipiens Becc, L. chinensis (Jacq.) R. Br. ex Mart., L. saribus (Lour.) Merr. ex A. Chev. (= L. cochinchinensis), L. subglobosa (Hassk.) Becc., Metroxylon sagu Rottb., Nipa sp., Oneosperma horrida, O. tigillarium (Jack) Ridl, Oreodoxa regia Kunth, Phoenix canariensis Chabaud, P. dactylifera L., P. sylvestris (L.) Roxb, P. theophrasti Greuter, Sabal umbraculifera (Jacq.) Mart.), Trachycarpus fortune (Hook.) H.Wendl. and Wash- ingtonia filifera Lindl.) H.Wendl.. Only one host plant species each was recorded from the 2 other families: Agave americana L. (Agavaceae) and Saccharum officinarum L. (Poaceae) (Murphy & Briscoe 1999; Faleiro 2006; Kontodimas et al. 2006; EPPO 2008; Malumphy & Molan 2009; Bertone et al. 2010; Dembilio & Jacas, 2012).

Fig. 3.

Close-up of model prediction across: A-Central and western United States and Mexico, B- Southeastern United States, C- North Africa, Southern Europe and Middle East and D- Caribbean region. White areas are deemed by the niche models to be unsuitable; gray areas and pink areas successively more suitable, and red areas to be quite suitable. Note the highly suitable areas in the region of Los Angeles in southern California, which coincide with the site where the species has successfully invaded.

Southern California as the newest colonization event of RPW and the Caribbean Islands were the impetus for this study. In California however, the suitable area is limited to a narrow fringe along the coast, with interior areas presenting only low suitabilities for the pest (Fig.3). As a consequence, our models suggest that, despite the current presence of the pest in Laguna Beach, California, the potential for direct spread in that state will be limited. Longer-distance transportation to other regions of the United States presenting suitable conditions (particularly in the Southeast) will need to be monitored carefully. Apart from California, no other reports exist for mainland South or North America; efforts should focus on preventing such establishment.

The spread of RPW across the world has accelerated since the middle 1980s. The original expectation regarding invasion pathways into North America was from the eastern side, perhaps coming from infested islands in the Caribbean. However, the first reported infestation was in California. Indeed, the morphology of the California invasive populations differs from that of specimens from Egypt, Europe, and the Caribbean (USDA 2010). Whether this difference reflects a separate invasion pathway, with independent origin from Asian populations, is unclear; however Hallet et al. (2004), based on morphological, molecular-genetic and breeding data considered R. ferrugineus and R. vulneratus (Panzer) as color morphs of the same species, and combined them under R. ferrugineus. Further studies need to link morphology (and potentially molecular characters as well) of the pest and its pathways for spread.

The cause of the high rate of spread of this pest has clearly been human intervention, by transporting infested young or adult date palm trees and offshoots from contaminated areas to uninfected areas (Abraham et al. 1998; Gomez & Ferry 2002). For instance, the introduction of the pest in the Caribbean in 2009 was the result of importation of date palms from Egypt to Aruba and Curaçao as part of a huge landscaping project undertaken by major tourist hotels. Similarly, introduction of the pest into Europe in 1993 occurred through importation of adult palms from Egypt to southern Spain (Gomez & Ferry 2002; Martin & Cabello 2006). In Egypt itself, introduction of the pest was via importation of offshoots from the United Arab Emirates (Ferry 1996). The knowledge of the situation in areas from where palms are purchased is therefore a very important step in reducing spread through trade and quarantine regulations should be enforced to ensure movement of RPW-free planting material.