Introduction
Species from one place can move to another due to anthropogenically driven activity and natural dispersal (Padayachee et al., Reference Padayachee, Irlich, Faulkner, Gaertner, Procheş, Wilson and Rouget2017). Alien species have the potential to become established and expand after being introduced to a new area, leading to ecological destruction and economic losses (Nahrung et al., Reference Nahrung, Liebhold, Brockerhoff and Rassati2023; Uden et al., Reference Uden, Mech, Havill, Schulz, Ayres, Herms, Hoover, Gandhi, Hufbauer, Liebhold and Marsico2023). There has been a dramatic rise in the prevalence of invasive species during the nineteenth century (Simberloff and Gibbons, Reference Simberloff and Gibbons2004). The invasiveness of these species is influenced by many factors, including climatic conditions, host plant resistance, propagule pressure and natural enemies in the invaded area (Skendžić et al., Reference Skendžić, Zovko, Pajač Živković, Lešić and Lemić2021).
Invasive pest population sizes, survival rates, geographic distributions, disease prevalence and proliferation can all be regulated by climatic factors (Aidoo et al., Reference Aidoo, Cunze, Guimapi, Arhin, Ablormeti, Tettey, Dampare, Afram, Bonsu, Obeng, Lutuf, Dickinson and Yankey2021, Reference Aidoo, Souza, da Silva, Santana, Picanço, Kyerematen, Sètamou, Ekesi and Borgemeister2022, Reference Aidoo, Souza, Silva, Júnior, Picanço, Heve, Duker, Ablormeti, Sétamou and Borgemeister2023a). Their populations, mainly in temperate regions, are expected to increase worldwide due to climate change (Schneider et al., Reference Schneider, Rebetez and Rasmann2022). However, it still needs to be determined how much of a role climate change plays in determining invasive species like the invasive Asian citrus psyllid Diaphorina citri Kuwayama (Hemiptera: Liviidae) in Ghana.
Diaphorina citri is a sap-sucking hemipteran pest of citrus and its close relatives. The psyllid undergoes five nymphal stages, which may take about 49 days, depending on the thermal conditions (Tsai and Liu, Reference Tsai and Liu2000). Apart from temperature, other factors that can influence the development and survival of the psyllid include availability of host plants and abundant flush growth. Diaphorina citri prefers warm and dry climates where conditions are suitable for developing all life stages (Hall et al., Reference Hall, Richardson, Ammar and Halbert2013).
The physical damage caused by the psyllid through feeding leads to the formation of sooty molds, and leaf distortion. Indirectly, the psyllid transmits the invasive pathogens; ‘Candidatus Liberibacter asiaticus’ (CLas) and ‘Ca. L. americanus’ (CLam) (Bové, Reference Bové2006). ‘Candidatus Liberibacter asiaticus’ occurs in Africa, Asia and the Americas, whereas CLam is distributed in Brazil. These phloem-limited pathogens have been associated with the deadliest disease of citrus called citrus greening disease or Huanglongbing (HLB). The interplay between the pathogens and the local accumulation of callose affects the movement of phloem in the plant (Welker et al., Reference Welker, Pierre, Santiago, Dutt, Vincent and Levy2022).
There is no effective method for preventing or curing HLB. Symptoms of HLB infection may not appear for months or even years, making it difficult to detect in affected trees (Lee et al., Reference Lee, Halbert, Dawson, Robertson, Keesling and Singer2015). Producing seedlings from CLas-free nursery stock, removing affected trees and applying insecticides to curb vectors are all part of HLB management. However, if pesticides are used extensively, pest resistance to pesticides may develop (Naeem et al., Reference Naeem, Freed, Jin, Akmal and Mehmood2016).
In 2023, for the first time, Ghana reported D. citri (Hemiptera: Liviidae) (Aidoo et al., Reference Aidoo, Ablormeti, Ninsin, Antwi-Agyakwa, Osei-Owusu, Heve, Dofuor, Soto, Edusei, Osabutey, Sossah, Aryee, Alabi and Sétamou2023b). The psyllids were obtained from orange jasmine Murraya paniculata (L.) in one of Ghana's 16 regions. However, none of the samples tested positive for the phloem-limited bacteria; CLas, CLam and CLaf. Risk maps of D. citri, which can serve as a guide for sampling and tracking, are urgently needed to develop quarantine and preventive measures. D. citri has been the subject of many species’ distribution models. For example, Aidoo et al. (Reference Aidoo, Souza, da Silva, Santana, Picanço, Kyerematen, Sètamou, Ekesi and Borgemeister2022) used the Maxent model to project where D. citri might develop in the world given the current and future climates. The potential distribution of D. citri was estimated in China using the Maxent and CLIMEX models (Wang et al., Reference Wang, Xiao and Zhang2015, Reference Wang, Yang, Wang, Zhang, Huang, Wen and Li2020). The Maxent model has a friendly user interface and works well when only presence and fewer occurrence records are available (Aidoo et al., Reference Aidoo, Tanga, Mohamed, Rasowo, Khamis, Rwomushana, Kimani, Agyakwa, Daisy, Sétamou, Ekesi and Borgemeister2019, Reference Aidoo, Souza, da Silva, Santana, Picanço, Kyerematen, Sètamou, Ekesi and Borgemeister2022).
The newly developed package maxnet in R is used to estimate Maxent models for regularised generalised linear models (Friedman et al., Reference Friedman, Hastie and Tibshirani2010). Maxnet employs the R package glmnet instead of relying on external sources such as the Maxent Java programme. It incorporates feature classes and regularization settings to effectively fit models comparable to those generated by the Maxent Java application (Phillips et al., Reference Phillips, Anderson, Dudík, Schapire and Blair2017). This enables enhanced integration of Maxent modelling with R's diverse range of visualization and analysis tools. The package helps to incorporate all the derived feature classes, with a particular emphasis on hinge features, as well as the default configured regularization settings from the Maxent Java programme into the maxnet package, thereby enabling the seamless and straightforward fitting of Maxent models within the R environment (Phillips et al., Reference Phillips, Anderson, Dudík, Schapire and Blair2017). The performance of maxnet is similar to that of the Maxent Java programmes (Phillips and Phillips, Reference Phillips and Phillips2021).
Understanding areas suitable for an invasive insect pest is a prerequisite for surveillance and tracking programmes. Using local distribution information of a species to examine its geographic distribution offers excellent accuracy when operating under the niche conservation tenet (Phillips, Reference Phillips2017). However, limiting occurrence records to only local records seems problematic given that D. citri has only recently become established in Ghana (Aidoo et al., Reference Aidoo, Ablormeti, Ninsin, Antwi-Agyakwa, Osei-Owusu, Heve, Dofuor, Soto, Edusei, Osabutey, Sossah, Aryee, Alabi and Sétamou2023b). These local records may not represent its ultimate niche in the country for accurate model predictions. Herein, we used the maxnet package in R to estimate the future spread risk for two shared socioeconomic pathways (SSPs) based on local and global occurrence records and climate datasets, and determine the climatic factors influencing the distribution of the pest. Also, we predicted the current suitable regions using local and global records to help policymakers develop effective responses to D. citri-borne disease.
Material and methods
Occurrence data
Four hundred and seventy (470) global occurrence records were sourced from Aidoo et al. (Reference Aidoo, Tanga, Mohamed, Rasowo, Khamis, Rwomushana, Kimani, Agyakwa, Daisy, Sétamou, Ekesi and Borgemeister2019), Aidoo et al. (Reference Aidoo, Ablormeti, Ninsin, Antwi-Agyakwa, Osei-Owusu, Heve, Dofuor, Soto, Edusei, Osabutey, Sossah, Aryee, Alabi and Sétamou2023b) and Sétamou et al. (Reference Sétamou, Soto, Tachin and Alabi2023). These records were supplemented by new field surveys in Ghana, and by data obtained from the Global Biodiversity Information Facility (GBIF). For the former, citrus and non-citrus host plants in the Western, Central, Eastern, Greater Accra, Ashanti and Volta Regions of Ghana were surveyed for the presence of the pest. If D. citri was found, Garmin eTrex® 32 × was used to obtain the Geographical Processing System (GPS) coordinates. The survey yielded 28 occurrence points. The occurrence records from GBIF were obtained using the {rgbif} package version 3.7.7 (Chamberlain et al., Reference Chamberlain, Barve, Mcglinn, Oldoni, Desmet, Geffert and Ram2023). Only records with a spatial resolution ≤1 km2 were retained for analysis. When only the names of locations were identified in the publications, their geographic coordinates were obtained online with Google Maps. The dataset was then formed by 498 occurrence points compiled from literature and field data, and 516 acquired from GBIF.
Occurrence records with a radius of 1 km around the capitals and centres of countries, equal absolute longitude and latitude, a radius of one degree around the GBIF headquarters, duplicated coordinates and zero values were removed, in addition to verifying that the informed coordinates corresponded to a valid coordinate reference system, using the clean_coordinate function of the {CoordinateCleaner} package version 2.0-20 (Zizka et al., Reference Zizka, Silvestro, Andermann, Azevedo, Duarte Ritter, Edler, Farooq, Herdean, Ariza, Scharn, Svantesson, Wengström, Zizka and Antonelli2019) (table S1). After these procedures, of the total 1014 points, 707 points of occurrence remained. After extracting environmental data in their coordinates and eliminating those in which at least one of the selected variables did not present information, the final set of occurrences was 704 points, of which 628 corresponded to invaded areas around the world, and others referred to D. citri native areas in South and Southeast Asia (e.g. Bangladesh, Bhutan, India, Nepal, Pakistan and southern China) (Bayles et al., Reference Bayles, Thomas, Simmons, Grafton-Cardwell and Daugherty2017; Bové, Reference Bové2014; Luo and Agnarsson, Reference Luo and Agnarsson2018).
To further mitigate the bias in occurrence data sampling, for the variable selection process, a spatial filter based on a 5 km distance was applied after data cleaning, considering the 2.5 km2 resolution of the used variables. For the final model, the 704 points of presence were filtered by applying an environmental filter (Varela et al., Reference Varela, Anderson, Garcia-Valdes and Fernandez-Gonzalez2014; Castellanos et al., Reference Castellanos, Huntley, Voelker and Lawing2019; Velazco et al., Reference Velazco, Svenning, Ribeiro and Laureto2020), using the function occfilt_env from the {flexsdm} package (Velazco et al., Reference Velazco, Rose, de Andrade, Minoli and Franklin2022). As the environmental filters are sensitive to the number of intervals (bins), filters with four, five, six, eight and ten bins were made, using the selected variables to build the distribution models. Spatial autocorrelation, based on the Moran Index, was calculated for each resulting dataset and the number of bins that showed the lowest autocorrelation (5 bins, with I = 0.352) was selected. After all these steps, a final set of 315 unique occurrence records was obtained, cleaned and filtered, for the modelling process (fig. 1).
Distribution of D. citri: points of occurrence (red, native areas; magenta, invaded areas) and study area (areas with darker contour).
As our intention was to estimate the areas susceptible to invasion by D. citri in Ghana, its probability of establishment was evaluated, and the results were interpreted without considering limitations to the future dispersal of the species, using native and invaded regions since this is the best option, especially if considering the modelling goal for invasive species (Beaumont et al., Reference Beaumont, Gallagher, Thuiller, Downey, Leishman and Hughes2009; Broennimann and Guisan, Reference Broennimann and Guisan2008; Zhang et al., Reference Zhang, Mammola, McLay, Capinha and Yokota2020).
Environmental data
The bioclimatic variables (table 1) utilised in this study were obtained from the Worldclim database version 2.1 (Fick and Hijmans, Reference Fick and Hijmans2017). These variables have an average spatial resolution of 2.5 arc-min, which corresponds to about 4 km at the equator. The data cover the period from 1970 to 2000. The {geodata} package version 0.5-8 (Hijmans et al., Reference Hijmans, Barbosa, Ghosh and Mandel2023) was used to access and process the data. These variables were chosen as predictor variables due to their ability to capture annual climate variations and limiting factors known to affect the geographic distribution of species (O'Donnell and Ignizio, Reference O'Donnell and Ignizio2012). In addition, Moderate Resolution Imaging Spectroradiometer (MODIS) Land Cover was sourced from Friedl and Sulla-Menashe (Reference Friedl and Sulla-Menashe2022).
Environmental variables used for the initial model
Variables that exhibited a Pearson's correlation coefficient (r) greater than the absolute value of 0.70, which is considered significant at an α level of 0.05, were categorised based on the results of a hierarchical cluster analysis (fig. 2). The cluster analysis was conducted using the {corrplot} package version 0.92 (Wei and Simko, Reference Wei and Simko2021). The initial set of all environmental variables was reduced using an algorithm that: (1) ranks variables based on permutation importance, (2) checks if the top-ranked variable is correlated with others (Pearson, r > |0.7|), and (3) if correlated, performs a Jackknife test between correlated variables, removing the variable that minimally impacts model performance when removed, based on the true skill statistics (TSS) metric (Allouche et al., Reference Allouche, Tsoar and Kadmon2006; Lawson et al., Reference Lawson, Hodgson, Wilson and Richards2014; Shabani et al., Reference Shabani, Kumar and Ahmadi2018). This process was repeated until all variables were tested, and all correlated variables were removed. Variables with low permutation importance (<3%) were then removed using Jackknife, ensuring that removal did not decrease model quality (based on estimated TSS value) using independent test data in ten permutations. We presented the contribution of the environmental variables to the base model in table S2. The resulting variables, and their descriptive statistics considering the points of presence, are presented in table 2. These variables were retained because findings from earlier studies indicate that, when developing models under various conditions, it is more useful to focus on a few variables that align with well-defined biological expectations rather than incorporating many variables with uncertain impacts on species distribution (Araújo and Guisan, Reference Araújo and Guisan2006; Austin and Van Niel, Reference Austin and Van Niel2011; Santini et al., Reference Santini, Benítez-López, Maiorano, Čengić and Huijbregts2021).
Correlation between bioclimatic variables. Blue colour sloping to the right indicates a positive correlation, while red sloping to the left indicates a negative correlation. The intensity of Pearson's correlation coefficient (r) increases from the circle (r = 0) to the ellipse (r = intermediate) to the line (r = 1). Correlated variables, r > 0.70, were grouped by centroid with hierarchical cluster analysis.
Descriptive statistics of the environmental variables considering their values at the 315 occurrence points (before the environmental filter) used to develop the models
Model development
All procedures related to data processing, development of models and maps were performed with the R environment, version 4.3.0 ‘Already Tomorrow’ (R Core Team, 2023) with the following packages: {rnaturalearth} version 1.0.1 (Massicotte and South, Reference Massicotte and South2023), to obtain spatial data from Ghana; {terra} version 1.7-29 (Hijmans, Reference Hijmans2023) and {sf} version 1.0-13 (Pebesma, Reference Pebesma2018), for analysis and transformation of spatial data; {flexsdm} version 1.3.3 (Velazco et al., Reference Velazco, Rose, de Andrade, Minoli and Franklin2022), for all species distribution modelling procedures, with resources from {maxnet} version 0.1.4 (Phillips and Phillips, Reference Phillips and Phillips2021); {pROC} version 1.18.2 (Robin et al., Reference Robin, Turck, Hainard, Tiberti, Lisacek, Sanchez and Müller2011), for graphics and estimates of the receiver operating characteristics (ROC) curve; {ade4} version 1.7-22 (Dray and Dufour, Reference Dray and Dufour2007) and {ecospat} version 3.5.1 (Broennimann et al., Reference Broennimann, Di Cola and Guisan2022), for ecological analyses; {factoextra} version 1.0.7 (Kassambara and Mundt, Reference Kassambara and Mundt2022), to visualise the results of multivariate analyses; {tmap} version 3.3-3 (Tennekes, Reference Tennekes2018), to plot all resulting maps.
The maximum entropy model (Maxent) was used because of its widespread application in species distribution modelling and its demonstrated effectiveness relative to alternative methods (Elith et al., Reference Elith, Graham, Anderson, Dudik, Ferrier, Guisan, Hijmans, Huettmann, Leathwick, Lehmann, Li, Lohmann, Loiselle, Manion, Moritz, Nakamura, Nakazawa, Overton, Peterson, Phillips, Richardson, Scachetti-Pereira, Schapire, Soberon, Williams, Wisz and Zimmermann2006, Reference Elith, Phillips, Hastie, Dudik, Chee and Yates2011; Heikkinen et al., Reference Heikkinen, Marmion and Luoto2012; Hijmans, Reference Hijmans2012; Venette, Reference Venette2017; Helmstetter et al., Reference Helmstetter, Conway, Stevens and Goldberg2021; Valavi et al., Reference Valavi, Guillera-Arroita, Lahoz-Monfort and Elith2022). Maxent has been employed to identify suitable regions for the establishment of D. citri on a global scale (Aidoo et al., Reference Aidoo, Souza, da Silva, Santana, Picanço, Kyerematen, Sètamou, Ekesi and Borgemeister2022). It has been recognised as a reliable approach for species distribution modelling (Valavi et al., Reference Valavi, Guillera-Arroita, Lahoz-Monfort and Elith2022). Yet, Maxent is susceptible to sample bias and can easily result in overfitting (Zhu et al., Reference Zhu, Liu and Gao2014). To mitigate the overfitting and enhance the generalizability of a model, it is crucial to carefully select parameters. This includes determining suitable transformations for the explanatory variables and identifying the best value for the regularization multiplier. This step was carried out through constructing and validating different models, seeking the best combination of these parameters.
The process of variable transformation in Maxent involves expanding the explanatory variables to include a broader range of derived variables. The bioclimatic variables chosen for the model are represented as derived variables, which are functions of the original variables (Phillips et al., Reference Phillips, Anderson and Schapire2006). This phenomenon might be interpreted as a change in the functional structure of the model specification, akin to including polynomial parts in a linear regression analysis. According to Phillips et al. (Reference Phillips, Anderson, Dudík, Schapire and Blair2017) and Phillips and Dudík (Reference Phillips and Dudík2008), the latest iterations of Maxent offer support for five types of transformations applicable to continuous independent variables. These transformations include linear, quadratic, threshold, forward hinge and reverse hinge.
A multi-step process (Perkins-Taylor and Frey, Reference Perkins-Taylor and Frey2020; Warren et al., Reference Warren, Matzke and Iglesias2020; Khan et al., Reference Khan, Li, Saqib, Khan, Habib, Khalid, Majeed and Tariq2022) was applied to identify species-specific adjustments used for optimizing the model (Shcheglovitova and Anderson, Reference Shcheglovitova and Anderson2013; Radosavljevic and Anderson, Reference Radosavljevic and Anderson2014). This was done to prevent excessive complexity and, consequently, low performance when projecting the model to different locations or under climate change scenarios (Elith et al., Reference Elith, Phillips, Hastie, Dudik, Chee and Yates2011; Merow et al., Reference Merow, Smith and Silander2013; Low et al., Reference Low, Zeng, Tan and Yeo2021). Several authors have employed a non-homogeneous Poisson process to estimate models (Renner and Warton, Reference Renner and Warton2013; Renner et al., Reference Renner, Elith, Baddeley, Fithian, Hastie, Phillips, Popovic and Warton2015; Phillips et al., Reference Phillips, Anderson, Dudík, Schapire and Blair2017). These models incorporated linear (L), quadratic (Q) and hinge (H) features and combinations of these features such as ‘LQ’, ‘QH’ and ‘LQH’. A base Maxent model was fitted using RM = 1 and FC = LQHP (default settings), employing the ‘Maxnet’ method and the training dataset (d), including all predictor variables. The goal was to evaluate and select the most important variables. A total of 119 models (FC = L, Q, H, LQ, QH, LQH, LQP, LQHP; RM = 1–5, with increments of 0.25) were fitted to identify the best combination of defined hyperparameters and to select the most suitable model (final model) using TSS with the selected variables (fig. S1). If results from different models are similar, the simplest one, ecologically easier to understand, was chosen. This is particularly important for the formulation of phytosanitary policies and decision-making processes.
Background points are based on occurrence distribution. In the first dataset we applied a geographic filter which is easy and faster, as we are working with all 19 variables. To run the final model, we filtered the occurrences with an environmental filter (just with the selected variables) and thus, had to create another background sample. For defining background points, the occurrence points were divided into four blocks using conventional k-fold cross-validation, allowing control over potential spatial autocorrelation between training and testing data and evaluating model transferability more appropriately (Roberts et al., Reference Roberts, Bahn, Ciuti, Boyce, Elith, Guillera-Arroita, Hauenstein, Lahoz-Monfort, Schröder, Thuiller, Warton, Wintle, Hartig and Dormann2017; Santini et al., Reference Santini, Benítez-López, Maiorano, Čengić and Huijbregts2021). Thirty grids with resolutions ranging from 0.5 to 8 degrees were generated, each with a minimum of ten occurrences (fig. S2 a, b). This method effectively addresses potential spatial autocorrelation between the training and test data, and provides a more suitable evaluation of the models’ transferability compared to alternative partitioning methods (Roberts et al., Reference Roberts, Bahn, Ciuti, Boyce, Elith, Guillera-Arroita, Hauenstein, Lahoz-Monfort, Schröder, Thuiller, Warton, Wintle, Hartig and Dormann2017; Santini et al., Reference Santini, Benítez-López, Maiorano, Čengić and Huijbregts2021). In their study on effects of future climate, land use and protected areas ineffectiveness on a dark scenario for Cerrado plant species, Velazco et al. (Reference Velazco, Villalobos, Galvão and De Marco Júnior2019) generated 20 grids, each with resolutions ranging from 0.5 to 5 degrees. The optimal grid size was selected considering (i) minimal spatial autocorrelation, (ii) maximum environmental similarity and (iii) minimal differences in records between training and testing data (Velazco et al., Reference Velazco, Rose, de Andrade, Minoli and Franklin2022). In this study, we made sure each cell included at least five occurrence points. For the purpose of conducting an autocorrelation test between groups, 50% of the points of existence were used. After establishing the division of occurrence records, the allocation of 10,000 background points was carried out in a random way, making sure they were dispersed proportionally based on the number of occurrences in each partition.
Excluding product features is justified by the fact that the model is defined by the marginal response curves of each predictor variable (Phillips et al., Reference Phillips, Anderson, Dudík, Schapire and Blair2017). These curves are more straightforward to interpret compared to those that rely on the values of other variables. It is worth noting that some of the used variables already encompass the combination of others. Phillips et al. (Reference Phillips, Anderson, Dudík, Schapire and Blair2017) specified the use of the clog-log output format as the most appropriate method for representing the probability of occurrence or establishment of a species.
The maximum entropy approach is a statistical technique that seeks to estimate an unknown probability distribution by leveraging the observed frequencies and background information within a defined study area. The goal of this strategy is to identify the distribution that shows the highest degree of geographic uniformity, often known as maximum entropy. This is achieved by incorporating the restrictions resulting from the environmental conditions seen in the known occurrence locations. Maxent often uses randomly generated background points distributed throughout the spatial extent of the study area (Phillips et al., Reference Phillips, Dudík, Elith, Graham, Lehmann, Leathwick and Ferrier2009). The purpose of this is to determine the extent of environmental variability surrounding the best conditions for the survival of a species. The number of background points significantly influences the model's predictions, and it is recommended to incorporate a considerable quantity of data to effectively represent the environmental conditions in which the species is found (Barbet-Massin et al., Reference Barbet-Massin, Jiguet, Albert and Thuiller2012; Northrup et al., Reference Northrup, Hooten, Anderson and Wittemyer2013). It is crucial to recognise that the results depict a model of habitat suitability relative in nature. Put otherwise, such findings indicate that one location is more suitable than another, without indicating the actual presence or inhabitation of the species.
The model development requires the spatial extent to encompass the accessible area for the species of interest during the relevant period (Araújo et al., Reference Araújo, Anderson, Márcia Barbosa, Beale, Dormann, Early, Garcia, Guisan, Maiorano, Naimi, O'Hara, Zimmermann and Rahbek2019). Additionally, the background data must be constrained to this same extent, and historical dispersion methods should be employed (Barve et al., Reference Barve, Barve, Jiménez-Valverde, Lira-Noriega, Maher, Peterson, Soberón and Villalobos2011; Merow et al., Reference Merow, Smith and Silander2013; Cooper and Soberón, Reference Cooper and Soberón2018). The choice of a proper geographic area for sampling background points varies depending on the species and the goals of the study (Santini et al., Reference Santini, Benítez-López, Maiorano, Čengić and Huijbregts2021). This is important for the development of niche models that rely on presence-only data, such as Maxent models (VanDerWal et al., Reference VanDerWal, Shoo, Johnson and Williams2009; Anderson and Raza, Reference Anderson and Raza2010; Barbet-Massin et al., Reference Barbet-Massin, Jiguet, Albert and Thuiller2012; Khosravi et al., Reference Khosravi, Hemani, Malekian, Flint and Flint2016; Cooper and Soberón, Reference Cooper and Soberón2018; Machado-Stredel et al., Reference Machado-Stredel, Cobos and Peterson2021; Amaro et al., Reference Amaro, Fidelis, da Silva and Marchioro2023). It is also crucial to make sure the sample size is large enough to adequately represent all environments (Renner et al., Reference Renner, Elith, Baddeley, Fithian, Hastie, Phillips, Popovic and Warton2015). In our study, the calibration area for the model was defined by considering the Köeppen–Geiger zones occupied by the species (Webber et al., Reference Webber, Yates, Le Maitre, Scott, Kriticos, Ota, McNeill, Le Roux and Midgley2011; Hill et al., Reference Hill, Gallardo and Terblanche2017), as depicted in fig. 3. The resulting calibration area encompassed 68,026,420 km2.
Map of the Köeppen–Geiger zones used to define the model calibration area, considering the current dispersion of Diaphorina citri.
Applying this approach is pertinent to models intended for extrapolation to different geographic regions beyond the calibration area or for alternative temporal intervals (Velazco et al., Reference Velazco, Rose, de Andrade, Minoli and Franklin2022). The best grid size was determined using the part_sblock function from the {flexsdm} package. This function systematically explores various block sizes and selects the most suitable size based on a multi-dimensional optimization procedure. The optimization process considers three dimensions: spatial autocorrelation (measured by Moran's I), environmental similarity (measured by Euclidean distance) and variation in data quantity among partition groups (measured by standard deviation – SD) (Velazco et al., Reference Velazco, Rose, de Andrade, Minoli and Franklin2022). Velazco et al. (Reference Velazco, Rose, de Andrade, Minoli and Franklin2022) established four partitions, each consisting of a test set of 40 cells and the resolving cells varying from 2 to 30 pixels.
In multiple cases of species distribution modelling, the chosen model predictions are transformed into binary maps that depict areas considered suitable or unsuitable for the species under investigation. While some scholars have advocated for the complete avoidance of this procedure unless there is a clear justification for applying the model (Guillera-Arroita et al., Reference Guillera-Arroita, Lahoz-Monfort, Elith, Gordon, Kujala, Lentini, McCarthy, Tingley and Wintle2015; Santini et al., Reference Santini, Benítez-López, Maiorano, Čengić and Huijbregts2021), the information provided by a map is valuable when addressing invasive species. It aids in focusing on areas for implementing phytosanitary public policies. The determination of a threshold value that serves as a cut-off point for site classification is a critical part, since varying thresholds might provide significantly divergent estimates (Jarnevich et al., Reference Jarnevich, Stohlgren, Kumar, Morisette and Holcombe2015). In this study, the threshold selected was max_sens_spec, which has been earlier shown to provide consistent findings (Liu et al., Reference Liu, Berry, Dawson and Pearson2005, Reference Liu, White and Newell2013, Reference Liu, Newell and White2016; Allouche et al., Reference Allouche, Tsoar and Kadmon2006). This threshold is equal to optimizing the vertical separation between a point on the ROC curve and the diagonal line, or maximizing the TSS. Still, it is crucial to acknowledge that when using thresholds to transform distribution models, the resulting interpretations should generally be regarded as hypothetical distributions (Liu et al., Reference Liu, White and Newell2013).
When constructing a model to estimate the possible distribution of a species, the model is adjusted to the observed data using a predetermined range of values for each environmental variable, known as the calibration range. When conducting model predictions for spatial or temporal projections, parts of the environment may show conditions that fall outside the calibration range. These conditions, under non-analogous circumstances, can lead to model extrapolation (Randin et al., Reference Randin, Dirnböck, Dullinger, Zimmermann, Zappa and Guisan2006; Williams et al., Reference Williams, Jackson and Kutzbach2007; Fitzpatrick and Hargrove, Reference Fitzpatrick and Hargrove2009; Owens et al., Reference Owens, Campbell, Dornak, Saupe, Barve, Soberón, Ingenloff, Lira-Noriega, Hensz, Myers and Peterson2013; Yates et al., Reference Yates, Bouchet, Caley, Mengersen, Randin, Parnell, Fielding, Bamford, Ban, Barbosa, Dormann, Elith, Embling, Ervin, Fisher, Gould, Graf, Gregr, Halpin, Heikkinen, Heinänen, Jones, Krishnakumar, Lauria, Lozano-Montes, Mannocci, Mellin, Mesgaran, Moreno-Amat, Mormede, Novaczek, Oppel, Ortuño Crespo, Peterson, Rapacciuolo, Roberts, Ross, Scales, Schoeman, Snelgrove, Sundblad, Thuiller, Torres, Verbruggen, Wang, Wenger, Whittingham, Zharikov, Zurell and Sequeira2018). To assess which areas had the highest extrapolation, the Shape metric (pers. comm.) was used, which consists of a model-independent method, since its calculation is based exclusively on the environmental distance between the training and projection data, that is, extrapolation is independent of model parameters and predictions.
Model evaluation
The scholarly literature has various perspectives on the best evaluation criteria for species distribution models. The ROC is a widely used metric for evaluating classification performance. One of its most commonly employed parts is the area under the curve (AUC) value. This metric is independent of threshold and scale (Merow et al., Reference Merow, Smith and Silander2013). Bradley (Reference Bradley1997) has revealed many useful attributes linked to the use of AUC as a performance metric in classification tasks. Still, the AUC values can be affected by the size of the study region (Smith, Reference Smith2013; Amaro et al., Reference Amaro, Fidelis, da Silva and Marchioro2023). Specifically, if the model encompasses a vast area and the species being studied has a limited distribution inside that area, the AUC values may be artificially inflated. Evaluating a classifier only based on total AUC may not be enough when assessing its performance specifically in areas characterised by high specificity or high sensitivity (Robin et al., Reference Robin, Turck, Hainard, Lisacek, Sanchez and Müller2009). To address these scenarios, the idea of partial area under the curve (pAUC) was created as a localised comparative method that specifically examines a segment of the ROC curve (McClish, Reference McClish1989; Jiang et al., Reference Jiang, Metz and Nishikawa1996; Streiner and Cairney, Reference Streiner and Cairney2007).
Alternative metrics for evaluating the adequacy of a model show certain limitations (Allouche et al., Reference Allouche, Tsoar and Kadmon2006). These limitations include sensitivity, which refers to the proportion of correctly predicted attendances or the true positive rate (TPR), specificity, which denotes the proportion of accurately predicted absences, and the TSS, a metric that remains unaffected by prevalence and is calculated using sensitivity and specificity (Allouche et al., Reference Allouche, Tsoar and Kadmon2006). Additionally, κ, another metric unaffected by prevalence, can also be used for this purpose (Allouche et al., Reference Allouche, Tsoar and Kadmon2006).
According to Hosmer et al. (Reference Hosmer, Lemeshow and Sturdivant2013), an ROC curve with an AUC value of 0.9 or higher indicates an excellent model fit. A value between 0.6 and 0.9 is considered good, while a value of 0.5 or lower suggests that the model performs no better than random chance. The TPR should be positioned near 1, hence indicating a heightened level of sensitivity. The TSS is a numerical measure that falls within the range of 0–1. A TSS value over 0.9 is regarded as ideal, while a value between 0.85 and 0.9 is considered exceptional. A TSS value ranging from 0.7 to 0.85 is classified as very good, while a value between 0.5 and 0.7 is considered good. Additionally, a TSS value ranging from 0.4 to 0.5 is considered decent, and any value ≤0.4 indicates a poor fit (Landis and Koch, Reference Landis and Koch1977). According to Peterson (Reference Peterson2006), to accurately represent low false-negative scores, omission values should be reduced and brought closer to zero.
Niche change
The occurrence records of the species were divided into two sets of data: the native area and the area invaded by D. citri. At first, updated Köppen–Geiger climate classes (Rubel and Kottek, Reference Rubel and Kottek2010) were extracted from a raster file obtained from the Climate Change and Infectious Diseases Group (https://koeppen-geiger.vu-wien.ac.at/present.htm) of the University of Vienna, Austria, to create histograms of the frequencies of occurrences and comparisons. To assess the ecological niches in the environment and compare the distinctions between the native and invaded regions, the COUE framework (Centroid shift, Overlap, Unfilling, and Expansion) proposed by Broennimann et al. (Reference Broennimann, Fitzpatrick, Pearman, Petitpierre, Pellissier, Yoccoz, Thuiller, Fortin, Randin, Zimmermann, Graham and Guisan2012) was used for both datasets. The comparison of niches was conducted in a two-dimensional space using the Schoener Index (D), a metric that measures the extent of overlap and ranges from 0 (showing no overlap) to +1 (representing complete overlap) (Broennimann et al., Reference Broennimann, Fitzpatrick, Pearman, Petitpierre, Pellissier, Yoccoz, Thuiller, Fortin, Randin, Zimmermann, Graham and Guisan2012). Subsequently, principal component analysis was used to analyse environmental data obtained from occurrence locations in both native and introduced areas to identify and assess the variability captured by the first two principal parts. The spatial extent of the environment was partitioned into a grid consisting of 100 cells, and the estimation of occurrence density within this spatial domain was conducted using a kernel density estimator, as described by Broennimann et al. (Reference Broennimann, Fitzpatrick, Pearman, Petitpierre, Pellissier, Yoccoz, Thuiller, Fortin, Randin, Zimmermann, Graham and Guisan2012) and Parravicini et al. (Reference Parravicini, Azzurro, Kulbicki and Belmaker2015).
The measured overlap (D) was assessed by using a niche equivalency test to determine whether it deviated significantly from a null distribution consisting of 1000 randomly generated D-metrics. This null distribution was created by randomly redistributing the occurrences of both niches among two datasets (Broennimann et al., Reference Broennimann, Fitzpatrick, Pearman, Petitpierre, Pellissier, Yoccoz, Thuiller, Fortin, Randin, Zimmermann, Graham and Guisan2012). The null hypothesis, which posits that the niches are comparable, is considered invalid when the observed D value falls below the fifth percentile of the null distribution. This null distribution was generated by considering the geographic availability of environmental conditions. Specifically, one niche was randomly distributed in the background while the other remained unchanged.
In this study, we used the same metrics as Guisan et al. (Reference Guisan, Petitpierre, Broennimann, Daehler and Kueffer2014) who examined various metrics related to niche dynamics in invasive species. These metrics include niche stability, which measures the proportion of environments within the introduced niche shared with the native niche. Another metric is niche expansion, which measures the proportion of environments in the introduced niche that do not overlap with the native niche. Last, the metric of unoccupied niche, also known as niche unfilling, assesses the proportion of environments that are not occupied by the invasive niche.
Projections for future climate scenarios
For future projections under different climate change scenarios, three periods were selected (2021–2040, 2041–2060 and 2061–2080), SSPs 245 and 585, using three global climate models (GCMs) considered suitable for Ghana (Oduro et al., Reference Oduro, Shuoben, Ayugi, Beibei, Babaousmail, Sarfo, Ullah and Ngoma2021), BCC-CSM2-MR, INM-CM5-0, MRI-ESM2-0, as per the Coupled Model Intercomparison Project Phase 6 (CMIP6).
The SSPs provide different development paths, contemplating possible trends about radiative forcing (W m–2). SSP 245 describes a society in which development follows a historical pattern without significant future deviations, with a radiative forcing in 2100 of 4.5 W m–2, representing an increase in global temperature between 1.4 and 2.8 °C. SSP 585 assumes a society in which economy is based on fossil fuels and intensive energy use, with a projected radiative forcing of 8.5 W m–2 in 2100 and a rise in global temperature between 3.5 and 5.5 °C.
The bioclimatic variables, SSPs and GCMs were obtained from the Worldclim. The mean model, for the periods and SSPs, was obtained by averaging the three GCMs used as a consensus model.
Contribution of environmental variables
Maxent's first formulation is equal to maximizing the probability of a parametric exponential distribution (Phillips et al., Reference Phillips, Dudík and Schapire2004), but more recently the identical maximum likelihood exponential model may be obtained through an Inhomogeneous Poisson Process (IPP) (Aarts et al., Reference Aarts, Fieberg and Matthiopoulos2012; Fithian and Hastie, Reference Fithian and Hastie2013; Renner and Warton, Reference Renner and Warton2013). One often used technique is the approach used by the ‘varImportance’ function within {fitMaxnet} package. The method used for calculating variable importance aligns with the method used in the R packages biomod2 and ecospat. The model generates predictions for every row in the combined environmental data table, created by stacking the occSWD and bkgSWD datasets. In the maxnet model object, the values of each variable undergo permutation across rows, resulting in a model prediction for each row using the permuted or shuffled data table. The permutation procedure is repeated numReplicates times for every variable. For every permutation, a Pearson correlation coefficient is calculated between the reference predictions and the predicted values obtained from the shuffled table. The significance value is calculated as the difference between 1-correlation coefficient. The outcome yields a vector of average scores for each variable, represented as a percentage relative to the sum of all average scores.
Results
Based on the TSS metric ( = 0.42224), the best combination was obtained using the linear, quadratic and hinge (LQH) classes simultaneously and regularised multiplication 0.5. The evaluation metrics of the selected model among the 50 tested models, provided by the {flexsdm} package and calculated from these, are presented in table 3, although all models have presented good performance to the random one.
Threshold-dependent and independent metrics (mean values and standard deviations) used to evaluate the developed model
The ROC curve of the final model (fig. 4a), resulting from evaluating true-positive predictions (sensitivity) and false-positive predictions (1 – specificity), showed a useful predictive capacity with an AUC value between 0.7 and 0.9. The ROC curve and partial AUC information when constraining the false- and true-positive rate (x = FPR = specificity; y = TPR = sensitivity) in the 90–100% range are illustrated in fig. 4b. The partial area (pAUC) can be interpreted as the mean sensitivity over the specified specificity range and the mean specificity over the specified sensitivity range. In addition, the suitability of habitats against the projection and training data are illustrated in fig. S3.
Graph of the receiver operating characteristic (ROC) curve for the potential geographic distribution model of Diaphorina citri in Ghana, showing the total area (AUC) under the ROC curve (a) and partial areas (pAUC) (b).
The response curves of the model are illustrated in fig. 5, which allow exploring the average marginal effect of environmental variables on the suitability of the environment for D. citri. The graphs show how the model response is individually influenced by each predictor variable, keeping the effects of the other variables fixed. The most suitable habitats (ideal values) for D. citri, as predicted by our model, are presented in supplemental information table S3. The histogram depicting the occurrence of D. citri against environmental variables is illustrated in fig. S4.
Graphs of the response curves of the variables used in the model. Mean diurnal range (Bio02), temperature seasonality (Bio04), mean temperature of warmest quarter (Bio10), precipitation seasonality (Bio15), precipitation of wettest quarter (Bio16), precipitation of driest quarter (Bio17), precipitation of coldest quarter (Bio19), Moderate Resolution Imaging Spectroradiometer (MODIS) Land Cover.
The analysis showed that the order of importance of the eight environmental variables is temperature seasonality (Bio04) > mean temperature of warmest quarter (Bio10) > precipitation of driest quarter (Bio17) > moderate resolution imaging spectroradiometer land cover (MODISLC) > precipitation seasonality (Bio15) > precipitation of coldest quarter (Bio19) > precipitation of wettest quarter (Bio16) > mean diurnal range (Bio02), were the most important drivers of D. citri distribution. Among these environmental variables, Bio04, Bio10, Bio17, MODISLC and Bio15 contributed to about 85% of the model (fig. 6). Further, the characteristics of the mean of the analysis on the D. citri are presented in table S4.
Percentage contribution of environmental variables to the base (a) and final (b) models. Annual mean temperature (Bio01), mean diurnal range (Bio02), isothermality (Bio03), temperature seasonality (Bio04), max temperature of warmest month (Bio05), min temperature of coldest month (Bio06), temperature annual range (Bio07), mean temperature of wettest quarter (Bio08), mean temperature of driest quarter (Bio09), mean temperature of warmest quarter (Bio10), mean temperature of coldest quarter (Bio11), annual precipitation (Bio12), precipitation of wettest month (Bio13), precipitation of driest month (Bio14), precipitation seasonality (Bio15), precipitation of wettest quarter (Bio16), precipitation of driest quarter (Bio17), precipitation of warmest quarter (Bio18), precipitation of coldest quarter (Bio19), Moderate Resolution Imaging Spectroradiometer (MODIS) Land Cover.
The potential geographic distribution of D. citri for Ghana, resulting from our model, is shown in fig. 7a, segmenting the probability of establishment into seven classes to help with visualization and comparison between locations. The high, optimal, moderate, marginal and unsuitable probability calculated at the occurrence points was 100% (table 4, fig. 7b). Applying a threshold that maximises the sum of sensitivity and specificity (max_sens_spec = 0.3656577) resulted in the map in fig. 7b which represents an area of 244,129 km2 (table 4). Extrapolation of models utilizing the shape metric is contingent upon its value, whereby a larger value (shown by shades of dark blue) indicates a greater disparity in the environment. Consequently, this discrepancy leads to a diminished level of reliability in the model's predictions for the corresponding area (fig. 7c).
The extrapolation map of the model, considering the presence and background points, elaborated from the extra_eval function of the {flexsdm} package, using the Shape metric. (A) Predicted probability of establishment of Diaphorina citri in Ghana. Warmer colours (closer to red) indicate environments more suitable for establishing the species. Small black dots refer to identified occurrences. (B) Regions considered optimal for the occurrence of D. citri in Ghana (area in red), considering max_sens_spec = 0.3656577, and (C) model extrapolation based on the Shape metric: the higher the value (shades of dark blue), the greater the environmental difference and, consequently, the lower the degree of reliability of a model prediction for the area.
Changes in habitat suitability from the current time to future with unsuitability, marginal, moderate, optimal and high percentages.
In its native distribution areas, D. citri predominantly occupies regions with hot semi-arid (steppe) climate (BSh), monsoon-influenced humid subtropical climate (Cwa) and tropical savanna, wet (Aw) climate classes according to the updated Köppen–Geiger climate classification (fig. 8a) with about 80% of the points belonging to these classes. Yet, in invaded areas, occurrences are more concentrated (about 58%) in the climate classes tropical savanna, wet (Aw), humid subtropical climate (Cfa) and tropical monsoon climate (Am) (fig. 8b). In general, we identified a greater number of climate classes occupied by D. citri in newly invaded compared to native areas, suggesting a change in niche.
Frequency histograms of the climatic classes occupied by the recorded occurrences of D. citri (native (A) and none-native (B) ranges, tropical monsoon climate (Am), tropical dry savanna climate (As), tropical savanna, wet (Aw), hot semi-arid (steppe) climate (Bsh), hot deserts climate (BWh), cold desert climate (BWk), humid subtropical climate (Cfa), monsoon-influenced humid subtropical climate (Cwa), subtropical highland climate or temperate oceanic climate with dry winters (Cwb), tropical rainforest climate (Af), cold semi-arid (steppe) climate (BSk), temperate oceanic climate (Cfb), hot-summer Mediterranean climate (Csa) and warm-summer Mediterranean climate (Csb).
The predicted distribution areas for D. citri in Ghana under the SSPs 245 and 585 and for the three periods are illustrated in fig. 9 and presented in table 4. The suitability classes for the country as predicted by the model are illustrated in fig. 10, with habitats ranging from marginal to high suitability. Yet, areas with high suitability are primarily located in the southern parts of the country. In citrus-producing regions, such as Ashanti, Central, Western, Eastern and Volta, a change in habitat suitability from the current time into the future is forecasted, with most parts showing high suitability for D. citri (fig. 10). The computed areas of habitat suitability for D. citri under the different climates are presented in table 4. We showed the optimal areas and the probability of establishment of the pest in the world in figs S5 and S6, respectively. Further, we illustrate Ghana's Maxent probability and classes maps in fig. S7. The probability of occurrence of D. citri across the time periods under climate change scenarios using a box plot are presented in fig. 11.
Mean future predictions of suitable areas for Diaphorina citri. (a) SSP 245 (2021-2040), (b) SSP 245 2041-2060, (c) SSP 245 2061-2080, (d) SSP 585 2021-2040, (e) SSP 585 2041-2060 and (f) SSP 585 2061-2080.
Classes of future predictions of suitable areas for Diaphorina citri. (a) SSP 245 (2021-2040), (b) SSP 245 2041-2060, (c) SSP 245 2061-2080, (d) SSP 585 2021-2040, (e) SSP 585 2041-2060 and (f) SSP 585 2061-2080.
Tukey climates from the current time (1970–2000) to the future (2021–2040; 2021–2060; 2061–2080) under SSPs 245 and 585.
Discussion
Linking climatic conditions to occurrence data is a commonly employed biogeographic strategy for characterizing species distributions and forecasting the potential effects of climate change (Guisan et al., Reference Guisan, Petitpierre, Broennimann, Daehler and Kueffer2014; Finch et al., Reference Finch, Butler, Runyon, Fettig, Kilkenny, Jose, Frankel, Cushman, Cobb, Dukes, Hicke and Amelon2021). Applying climate modelling in assessing habitat suitability has clearly demonstrated that climate change will substantially affect the distribution patterns of invasive species (Aidoo et al., Reference Aidoo, Souza, da Silva, Santana, Picanço, Kyerematen, Sètamou, Ekesi and Borgemeister2022, Reference Aidoo, Souza, Silva, Júnior, Picanço, Heve, Duker, Ablormeti, Sétamou and Borgemeister2023a). Notwithstanding the limitations and uncertainties in the outcomes of species distribution models, using these models is a valuable technique for predicting potential distribution of species (Woodman et al., Reference Woodman, Forney, Becker, DeAngelis, Hazen, Palacios and Redfern2019). Previous studies have established a correlation between temperature and precipitation patterns and their impacts on the spread, survival and development of D. citri and other invasive species (e.g. Hall et al., Reference Hall, Hentz, Meyer, Kriss, Gottwald and Boucias2012; Devi and Sharma, Reference Devi and Sharma2014). The present work used the CMIP6 data to predict the potential geographic distribution of D. citri in Ghana. Currently, the pest has a broad geographical distribution within the country, encompassing the major citrus-growing regions including, Western, Eastern, Central, Volta, Northern and the Ashanti Region (Brentu et al., Reference Brentu, Oduro, Offei, Odamtten, Vicent, Peres and Timmer2012; Asare-Bediako et al., Reference Asare-Bediako, Addo-Quaye, Tetteh, Buah, Van Der Puije and Acheampong2013). Our model now predicts a further expansion, with new areas of habitat suitability covering parts of Ahafo, Bono, Bono East, Northern, Savannah, North East, Upper East, Oti and Western North Regions. However, under the current climate areas with high and optimal suitability are restricted to southern parts of the country, including parts of Volta, Ashanti, Eastern, Western, Western North, Greater Accra and Central Regions, with moderate to marginal suitable climates in parts of the more arid Northern, North East, Upper West and Upper East Regions of Ghana. Our modelling results suggest that about 244,129 km2 of Ghana are climatically suitable for D. citri, which corresponds well with the available occurrence data (Aidoo et al., Reference Aidoo, Ablormeti, Ninsin, Antwi-Agyakwa, Osei-Owusu, Heve, Dofuor, Soto, Edusei, Osabutey, Sossah, Aryee, Alabi and Sétamou2023b), indicating that the maxnet package can effectively predict the pest's habitats. The original formulation of Maxent equals maximizing the likelihood of a parametric exponential distribution (Phillips et al., Reference Phillips, Dudík and Schapire2004). Like the Maximum Entropy java software, the maxnet package in R can generate forecasts and decrease commission mistakes when only a few occurrence records for a species are accessible.
In this study, temperature seasonality, mean temperature of warmest quarter, precipitation of driest quarter, moderate resolution imaging spectroradiometer land cover, precipitation seasonality, precipitation of coldest quarter, precipitation of wettest quarter and mean diurnal range were identified as significant environmental factors, collectively accounting for more than two-thirds of the observed geographical distribution of the pest. Earlier research on the growth and spread of D. citri has shown temperature and precipitation are crucial in influencing its reproductive, developmental, migratory, morphology and dispersal patterns (Antolínez et al., Reference Antolínez, Olarte-Castillo, Martini and Rivera2022; Paris et al., Reference Paris, Allan, Hall, Hentz, Croxton, Ainpudi and Stansly2017). Further, temperature alters the flight capacity of D. citri (Antolínez et al., Reference Antolínez, Olarte-Castillo, Martini and Rivera2022). New citrus flush production as well as maximum temperature, daily lowest temperature and rainfall have been positively correlated with D. citri infestations (Zorzenon et al., Reference Zorzenon, Tomaseto, Daugherty, Lopes and Miranda2021). The same authors found that the migration patterns of D. citri aligned with seasonal variations in specific climatic factors. Higher levels of humidity and daily maximum temperatures were associated with adverse effects, while increased rainfall in the preceding weeks had a favourable impact. Specifically, ideal conditions for the psyllid's spread include a mean diurnal range ranging from 2.7 to 4.2 °C, temperature seasonality of 56.5 °C and a mean temperature of warmest quarter of 22.2 °C (table S3). The changes in the climatic condition pattern resulting from climate change have the potential to exert both direct and indirect effects on the survival of invasive species. Nevertheless, it should be noted that in Brazil, the optimal environmental parameters for D. citri are typically observed within a temperature range of 22–28 °C (Zavala-Zapata et al., Reference Zavala-Zapata, Lázaro-Dzul, Sánchez-Borja, Vargas-Tovar, Álvarez-Ramos and Azuara-Domínguez2022). However, laboratory studies predictions showed that temperatures ranging from 25 to 28 °C were optimal for the proliferation of D. citri and the subsequent expansion of its population (Tsai and Liu, Reference Tsai and Liu2000). Previous research conducted under controlled and standardised thermal conditions demonstrated that D. citri is unable to complete its life cycle above a temperature barrier of 33 °C (Tsai and Liu, Reference Tsai and Liu2000). Yet, a more recent investigation by Milosavljević et al. (Reference Milosavljević, McCalla, Morgan and Hoddle2020) has provided empirical evidence that D. ciri has the ability to undergo developmental processes at temperatures exceeding 35 °C, under both constant and fluctuating temperature regimes. Furthermore, Antolínez et al. (Reference Antolínez, Olarte-Castillo, Martini and Rivera2022) demonstrated an effective completion of the psyllid's life cycle throughout a range of temperature treatments spanning from 28 to 43 °C. However, it is important to note that the daily cycles with a temperature of 43 °C sustained for 6 h posed an exception to this pattern. Moreover, these authors observed a significant decline in the rate of adult emergence as the temperature surpassed 40 °C. In their study, Hall et al. (Reference Hall, Wenninger and Hentz2011) conducted an estimation of temperature thresholds pertaining to the oviposition behaviour of D. citri. The results indicated that the lower and higher thresholds for oviposition were 16 and 41.6 °C, respectively. Comparable ranges for D. citri response curves were found for the warmest quarter's mean temperature, warmest month's maximum temperature and temperature seasonality: 24.46–34.27, 28.6–40.91 and 56.83–818.03 °C, respectively (Wang et al., Reference Wang, Yang, Luo, Wang, Lu, Huang, Zhao and Li2019). According to Skendžić et al. (Reference Skendžić, Zovko, Pajač Živković, Lešić and Lemić2021), global climate warming has the potential to affect insects, including the expansion of their geographic ranges, enhanced survival during winter months, an increase in the number of generations, heightened susceptibility to invasive insect species and insect-borne plant diseases and alterations in the dynamics of their interactions with host plants and natural predators.
Our predictions showed that precipitation seasonality of 121.7 mm, precipitation of wettest quarter of 25.9–317.1 mm, precipitation of driest quarter of 459.5 mm and precipitation of coldest quarter of 183.8 mm were ideal for the pest's proliferation in Ghana. Smaller insects such as D. citri are very vulnerable to severe precipitation which can dislodge or wash them off from their hosts. For example, Beattie (Reference Beattie2020) showed that heavy rains washed off eggs and nymphs, thereby considerably affecting D. citri populations. The forecast conducted in China indicated that the optimal range of precipitation during the wettest quarter for D. citri was found to be between 562.89 and 1189.75 mm. Similarly, the range of precipitation during the warmest quarter was seen to vary from 503.73 to 1533.58 mm (Wang et al., Reference Wang, Yang, Luo, Wang, Lu, Huang, Zhao and Li2019). Several studies have shown that the primary constraint on the population size, geographical range and possible spread of D. citri is the cold temperatures experienced during winter months (e.g. Hall et al., Reference Hall, Wenninger and Hentz2011; López-Collado et al., Reference López-Collado, López-Arroyo, Robles-García and Márquez-Santos2013). This suggests that the survival and distribution of D. citri are primarily influenced by climatic factors such as temperature and precipitation.
Results from the maxnet package indicated high-suitability habitats for D. citri, which overlap with the major citrus-producing regions in Ghana. Hence, the model showed a high level of reliability in forecasting the distribution of D. citri confirming previous findings (Aidoo et al., Reference Aidoo, Ablormeti, Ninsin, Antwi-Agyakwa, Osei-Owusu, Heve, Dofuor, Soto, Edusei, Osabutey, Sossah, Aryee, Alabi and Sétamou2023b). Based on the model outputs for the projected future climatic conditions of the SSPs245 and 585 from the 2040s to the 2080s, Ghana's regions including Volta, Greater Accra, Central, Western and Eastern constitute areas that harbour highly suitable climatic conditions for D. citri. In contrast, some of Ghana's northern regions, including Savannah, Northern, North East, Upper East and Upper West, are projected to be less suitable for the psyllid. The climate suitable areas will increase from the current time until the 2080s for both climate change scenarios, with a few areas in the future showing unsuitability for the pest.
A study by Aidoo et al. (Reference Aidoo, Ablormeti, Ninsin, Antwi-Agyakwa, Osei-Owusu, Heve, Dofuor, Soto, Edusei, Osabutey, Sossah, Aryee, Alabi and Sétamou2023b) reported for the first time the presence of D. citri in the Volta Region of Ghana and suggested regular tracking and surveillance of the pest. Our risk maps can serve as a guide for the future development of control measures to successfully manage the pest under current and future climate change conditions. Such pest management measures are urgently needed to prevent the further spread of D. citri in Ghana and beyond because of the expected rising temperatures expanding suitability for the pest in most parts of the country. However, even a reduced distribution of D. citri could still threaten citrus production in Ghana. Moreover, invasive species can often adapt to changing climates (Barrett, Reference Barrett, Mooney and Hobbs2000). Hence, the national phytosanitary authorities in Ghana should continue their efforts in preventing the spread of D. citri in the country as well as potential spillovers into the region.
The response curves from our model illustrate the manner in which the predicted probability of the species’ habitat alters in response to fluctuations in each predictor, while holding all other variables at their respective average sample values. They offer insights into the environmental factors correlated with varying species occurrence probabilities (Merow et al., Reference Merow, Smith and Silander2013; Tesfamariam et al., Reference Tesfamariam, Gessesse and Melgani2022). These parameters may correspond to the species’ established ecological preferences or tolerances, explaining certain features of its biology. However, species’ biology is frequently shaped by many interconnected elements, encompassing biotic interactions, reproductive strategies and other relevant characteristics. The ecological significance of response curves may vary depending on the particular context and size of the investigation. Varying reaction patterns may arise due to different life phases, geographical regions or year's seasons. Maxent models can capture intricate interactions; however, it is essential to note that overfitting, which entails fitting noise present in the data, can lead to unrealistic response curves that then need to be revised (Phillips et al., Reference Phillips, Dudík, Elith, Graham, Lehmann, Leathwick and Ferrier2009). Thorough model validation and careful evaluation of ecological plausibility are crucial aspects to be considered when modelling species distribution. Nonetheless, to thoroughly comprehend species biology, it is necessary to incorporate broader ecological information, field observations and domain experience alongside the findings obtained using maxnet in future investigations and analysis.
SDMs, like the maxnet package in R, can have some important drawbacks. We did not consider the potential effects of socioeconomic development, D. citri evolution and adaptation, the introduction of new plant protection and regulatory services policies, and changed farm-level management practices, all of which could considerably change the spread risk, in addition to human movement and interventions. Our predictions are based on SSPs, and the chosen combinations of emission and socioeconomic scenarios can influence our projections. Thus, different modelling methods and their combinations should be used in future studies to provide a better understanding of the uncertainty around our estimations. Biotic and abiotic factors such as pest pressure, predators, parasitoids and elevation that were not included in our model should be considered in future predictions. Notwithstanding these drawbacks, our results offer valuable data for developing surveillance and preventive policies against a further spread of D. citri in Ghana and the pest's impact on the country's citrus sector.
Conclusion
The present investigation examined the suitable ecological environment for D. citri in Ghana, including both present and projected climate change conditions represented by the SSP245 and SSP585 scenarios. Such an analysis is crucial for the formulation of effective strategies and policies aimed at an effective control of D. citri. The potential habitats of the pest was predicted using the maxnet in R. Our findings indicate that temperature and rainfall conditions are significant predictors for the possible distribution of D. citri in Ghana. Thus, the current climatic regions in Ghana are highly conducive for D. citri and are expected to expand until the 2080s. Still, it is projected that certain crucial citrus-producing regions in southern Ghana will continue to exhibit a high level of suitability for D. citri. Our findings can thus assist researchers and policymakers in developing effective and well-targeted pest management strategies for D. citri in a changing climate.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0007485324000105.
Data availability
The datasets generated in this study are available from the corresponding author upon request.
Acknowledgements
We thank the Universidade Federal dos Vales do Jequitinhonha e Mucuri (UFVJM) and the University of Environment and Sustainable Development (UESD) for access to instruments. We would like to thank the National Council for Scientific and Technological Development (Conselho Nacional de Desenvolvimento Científico e Tecnológico – CNPq) and the Minas Gerais State Research Foundation (Fundação de Amparo a Pesquisa do Estado de Minas Gerais – FAPEMIG) and the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior – Brasil (CAPES) – Finance Code 001.
Author contributions
The manuscript was conceived, designed and written by K. D. N. and O. F. A. Material preparation and data collection were conducted by K. D. N., G. C. A., J. O.-O., O. F. A., F. K. A., A. K. D., W. K. H., G. E., L. K. A., P. B. and H. A. B. Analysis were performed by P. G. C. S., G. C. A., E. J. D. V. B. and O. F. A. The manuscript was reviewed and edited by R. S. d. S. C. B. and M. S. contributed to review and editing of the article. All authors read and approved the final manuscript.
Financial support
Universidade Federal dos Vales do Jequitinhonha e Mucuri (UFVJM).
Competing interests
None.
