Risk of spread of the Asian citrus psyllid Diaphorina citri Kuwayama (Hemiptera: Liviidae) in Ghana

Abstract

The impact of invasive species on biodiversity, food security and economy is increasingly noticeable in various regions of the globe as a consequence of climate change. Yet, there is limited research on how climate change affects the distribution of the invasive Asian citrus psyllid Diaphorina citri Kuwayama (Hemiptera:Liviidae) in Ghana. Using maxnet package to fit the Maxent model in R software, we answered the following questions; (i) what are the main drivers for D. citri distribution, (ii) what are the D. citri-specific habitat requirements and (iii) how well do the risk maps fit with what we know to be correctly based on the available evidence?. We found that temperature seasonality (Bio04), mean temperature of warmest quarter (Bio10), precipitation of driest quarter (Bio17), moderate resolution imaging spectroradiometer land cover and precipitation seasonality (Bio15), were the most important drivers of D. citri distribution. The results follow the known distribution records of the pest with potential expansion of habitat suitability in the future. Because many invasive species, including D. citri, can adapt to the changing climates, our findings can serve as a guide for surveillance, tracking and prevention of D. citri spread in Ghana.

Introduction

Species from one place can move to another due to anthropogenically driven activity and natural dispersal (Padayachee et al., 2017). 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., 2023; Uden et al., 2023). There has been a dramatic rise in the prevalence of invasive species during the nineteenth century (Simberloff and Gibbons, 2004). 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., 2021).

Invasive pest population sizes, survival rates, geographic distributions, disease prevalence and proliferation can all be regulated by climatic factors (Aidoo et al., 2021, 2022, 2023a). Their populations, mainly in temperate regions, are expected to increase worldwide due to climate change (Schneider et al., 2022). 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, 2000). 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., 2013).

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é, 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., 2022).

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., 2015). 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., 2016).

In 2023, for the first time, Ghana reported D. citri (Hemiptera: Liviidae) (Aidoo et al., 2023b). 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. (2022) 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., 2015, 2020). The Maxent model has a friendly user interface and works well when only presence and fewer occurrence records are available (Aidoo et al., 2019, 2022).

The newly developed package maxnet in R is used to estimate Maxent models for regularised generalised linear models (Friedman et al., 2010). 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., 2017). 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., 2017). The performance of maxnet is similar to that of the Maxent Java programmes (Phillips and Phillips, 2021).

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, 2017). However, limiting occurrence records to only local records seems problematic given that D. citri has only recently become established in Ghana (Aidoo et al., 2023b). 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. (2019), Aidoo et al. (2023b) and Sétamou et al. (2023). 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., 2023). Only records with a spatial resolution ≤1 km² 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., 2019) (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., 2017; Bové, 2014; Luo and Agnarsson, 2018).

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 km² resolution of the used variables. For the final model, the 704 points of presence were filtered by applying an environmental filter (Varela et al., 2014; Castellanos et al., 2019; Velazco et al., 2020), using the function occfilt_env from the {flexsdm} package (Velazco et al., 2022). 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).

Figure 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., 2009; Broennimann and Guisan, 2008; Zhang et al., 2020).

Environmental data

The bioclimatic variables (table 1) utilised in this study were obtained from the Worldclim database version 2.1 (Fick and Hijmans, 2017). 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., 2023) 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, 2012). In addition, Moderate Resolution Imaging Spectroradiometer (MODIS) Land Cover was sourced from Friedl and Sulla-Menashe (2022).

Table 1. 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, 2021). 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., 2006; Lawson et al., 2014; Shabani et al., 2018). 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, 2006; Austin and Van Niel, 2011; Santini et al., 2021).

Figure 2. 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.

Table 2. 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, 2023), to obtain spatial data from Ghana; {terra} version 1.7-29 (Hijmans, 2023) and {sf} version 1.0-13 (Pebesma, 2018), for analysis and transformation of spatial data; {flexsdm} version 1.3.3 (Velazco et al., 2022), for all species distribution modelling procedures, with resources from {maxnet} version 0.1.4 (Phillips and Phillips, 2021); {pROC} version 1.18.2 (Robin et al., 2011), for graphics and estimates of the receiver operating characteristics (ROC) curve; {ade4} version 1.7-22 (Dray and Dufour, 2007) and {ecospat} version 3.5.1 (Broennimann et al., 2022), for ecological analyses; {factoextra} version 1.0.7 (Kassambara and Mundt, 2022), to visualise the results of multivariate analyses; {tmap} version 3.3-3 (Tennekes, 2018), 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., 2006, 2011; Heikkinen et al., 2012; Hijmans, 2012; Venette, 2017; Helmstetter et al., 2021; Valavi et al., 2022). Maxent has been employed to identify suitable regions for the establishment of D. citri on a global scale (Aidoo et al., 2022). It has been recognised as a reliable approach for species distribution modelling (Valavi et al., 2022). Yet, Maxent is susceptible to sample bias and can easily result in overfitting (Zhu et al., 2014). 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., 2006). 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. (2017) and Phillips and Dudík (2008), 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, 2020; Warren et al., 2020; Khan et al., 2022) was applied to identify species-specific adjustments used for optimizing the model (Shcheglovitova and Anderson, 2013; Radosavljevic and Anderson, 2014). 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., 2011; Merow et al., 2013; Low et al., 2021). Several authors have employed a non-homogeneous Poisson process to estimate models (Renner and Warton, 2013; Renner et al., 2015; Phillips et al., 2017). 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., 2017; Santini et al., 2021). 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., 2017; Santini et al., 2021). 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. (2019) 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., 2022). 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., 2017). 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. (2017) 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., 2009). 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., 2012; Northrup et al., 2013). 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., 2019). Additionally, the background data must be constrained to this same extent, and historical dispersion methods should be employed (Barve et al., 2011; Merow et al., 2013; Cooper and Soberón, 2018). 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., 2021). This is important for the development of niche models that rely on presence-only data, such as Maxent models (VanDerWal et al., 2009; Anderson and Raza, 2010; Barbet-Massin et al., 2012; Khosravi et al., 2016; Cooper and Soberón, 2018; Machado-Stredel et al., 2021; Amaro et al., 2023). It is also crucial to make sure the sample size is large enough to adequately represent all environments (Renner et al., 2015). 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., 2011; Hill et al., 2017), as depicted in fig. 3. The resulting calibration area encompassed 68,026,420 km².

Figure 3. 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., 2022). 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., 2022). Velazco et al. (2022) 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., 2015; Santini et al., 2021), 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., 2015). In this study, the threshold selected was max_sens_spec, which has been earlier shown to provide consistent findings (Liu et al., 2005, 2013, 2016; Allouche et al., 2006). 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., 2013).

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., 2006; Williams et al., 2007; Fitzpatrick and Hargrove, 2009; Owens et al., 2013; Yates et al., 2018). 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., 2013). Bradley (1997) 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, 2013; Amaro et al., 2023). 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., 2009). 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, 1989; Jiang et al., 1996; Streiner and Cairney, 2007).

Alternative metrics for evaluating the adequacy of a model show certain limitations (Allouche et al., 2006). 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., 2006). Additionally, κ, another metric unaffected by prevalence, can also be used for this purpose (Allouche et al., 2006).

According to Hosmer et al. (2013), 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, 1977). According to Peterson (2006), 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, 2010) 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. (2012) 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., 2012). 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. (2012) and Parravicini et al. (2015).

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., 2012). 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. (2014) 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., 2021), 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., 2004), but more recently the identical maximum likelihood exponential model may be obtained through an Inhomogeneous Poisson Process (IPP) (Aarts et al., 2012; Fithian and Hastie, 2013; Renner and Warton, 2013). 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.

Table 3. 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.

Figure 4. 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.

Figure 5. 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.

Figure 6. 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 km² (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).

Figure 7. 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.

Table 4. 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.

Figure 8. 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.

Figure 9. 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.

Figure 10. 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.

Figure 11. 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., 2014; Finch et al., 2021). 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., 2022, 2023a). 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., 2019). 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., 2012; Devi and Sharma, 2014). 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., 2012; Asare-Bediako et al., 2013). 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 km² of Ghana are climatically suitable for D. citri, which corresponds well with the available occurrence data (Aidoo et al., 2023b), 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., 2004). 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., 2022; Paris et al., 2017). Further, temperature alters the flight capacity of D. citri (Antolínez et al., 2022). 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., 2021). 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., 2022). 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, 2000). 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, 2000). Yet, a more recent investigation by Milosavljević et al. (2020) 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. (2022) 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. (2011) 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., 2019). According to Skendžić et al. (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 (2020) 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., 2019). 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., 2011; López-Collado et al., 2013). 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., 2023b). 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. (2023b) 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, 2000). 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., 2013; Tesfamariam et al., 2022). 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., 2009). 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.