Seasonal forecasting of Bactrocera dorsalis Hendel, 1912 (Diptera: Tephritidae) in bioclimatic zones of Sri Lanka using the SARIMA model

Bactrocera dorsalis Hendel is a severe fruit pest that causes significant economic losses globally. Despite B. dorsa-lis having been distributed mostly across Asia, studies on its current and future density variation in Sri Lanka are sparse to date. The present study was thus carried out to assess the contemporary density variation (2020–2022) and future density fluctuation (2023–2025) of B. dorsalis in bioclimatic zones of Sri Lanka. The density was assessed using the monthly-based fruit fly trap collection method from randomly selected 40 locations in all bioclimatic zones (wet, intermediate, dry, and arid). The SARIMA modelling technique was applied for delineating the best-fit model and for density forecasting in each bioclimatic zone. The density variations were depicted for the year and for the bioclimatic zone (2020–2025) by colour intensity maps using QGIS. According to the findings, B. dorsalis shows a seasonal component to its year-round density variation and an ascending trend in its density from 2020 to 2025. Density forecasting records a 20%, 30%, 26%, and 37% density increase in the wet, intermediate, dry, and arid zones, respectively, in 2025. In 2025, the highest predicted B. dorsalis density from the arid zone and the lowest predicted density from the wet zone were recorded. This study contains the first forecasting attempt for B. dorsalis density using the SARIMA approach as well as the application of colour-intensity depiction for its density variation in Sri Lanka, which leads decision makers and stakeholders in economic agriculture to plan the scientific management of B. dorsa-lis to avoid its current and potential future threat to the country’s fruit industry.


Introduction
Fruit flies (Diptera: Tephritidae) are a group of important insect pest.Most fruit flies are polyphagous, and they cause severe economic loss by damaging many horticulture products, especially fruits and vegetables, in the tropical region (White and Elson 1994;Clarke et al. 2005;Wei et al. 2019;Peng et al. 2020).Among them, B. dorsalis is considered a severe fruit pest worldwide (Wei et al. 2019;Peng et al. 2020), as well as being classified as a level 1 quarantine pest by many countries (Liu et al. 2019;CABI 2021).This species has a very high invasive capability, and hence its damage has been reported up to 100% for many economically important fruit types, resulting in huge economic loss (Ekesie et al. 2009;Liendo et al. 2018;CABI 2021).A widespread distribution of B. dorsalis is reported throughout Asia, with limited distribution in Pacific and African regions (Aketarawong et al. 2007;Drew and Romig 2013;Choudhary et al. 2016).
Sri Lanka is a tropical country, and its climatic conditions are favourable for high floral and faunal diversity.Based on the climatic factors, Sri Lanka is divided into four main bioclimatic zones: wet, intermediate, dry, and arid.The distribution of B. dorsalis as the predominant fruit fly species in all bioclimatic zones of Sri Lanka was reported by Wijekoon et al. (2023).
In several countries outside of Sri Lanka, forecasting models have been extensively used for the prediction of tephritid populations (Worner 1988;Yonow and Sutherst 1998;Sutherst et al. 2000;Vera et al. 2002;Stephens et al. 2007), and especially for B. dorsalis density by Dong et al. (2022).SARIMA (Seasonal Autoregressive Moving Average) is a seasonal time series forecasting model that is commonly used for studying spatial variability through time (Andrew 1994).This model has been extensively applied to forecast population growth with a seasonal component.The Box-Jenkins methodology is referred to as the systematic method to identify, fit, and check the data set in the SARIMA time series models (Cooray 2008).This model is scientifically vital to predict insect pest density, which is of agro-ecological significance.Though the application of SARIMA forecasting of time series has been reported in several studies in other countries (Saboia 1997;He and Lv 1992;Fan and Lv 1997;Meng 1997;Yang 2008;Zeng and Quanpeng 2011;Sanchez et al. 2013;Adebiyi et al. 2014;Unnikrishnan and Suresh 2016;Wang et al. 2016;Guo et al. 2020;Sun et al. 2020;Duan 2021;Hu & Wu 2021;Lin-feng 2021), few records have been reported in Asia (Narava et al. 2022).
To date, no study has been carried out to determine the density variation and forecasting of B. dorsalis in Sri Lanka.As the first attempt in Sri Lanka, the use of the SARIMA model to forecast B. dorsalis population is vital for agriculture stakeholders to determine the advent and patterns of B. dorsalis density over time and the timely intervention of proper pest control methods before damage occurs.Moreover, the outputs of the SARIMA application will be helpful to comprehend the contest of B. dorsalis behaviour and its potential future outbreak patterns in spatial scenarios.The present study was thus intended to assess the seasonal variation, modelling, and future forecasting of B. dorsalis density in bioclimatic zones of Sri Lanka.

Study area
Sri Lanka is an island with an area of 65,610 km 2 and is located in the tropical region.This study was conducted in all four bioclimatic zones (wet, intermediate, dry, and arid) of Sri Lanka and these four major bioclimatic zones have been divided based on rainfall, temperature, and relative humidity.

Site selection and sample collection
Ten sampling locations in each bioclimatic zone and 40 locations from all four bioclimatic zones were randomly selected.The GIS (Geographical Information System) coordination data for each study location was recorded.The selected study locations are depicted in Fig. 1, and their GIS coordination is addressed in Table 1.In each location, 100 m 2 of area were selected, and a ME (methyleugenol) (5 cm diameter, 10 cm height, two circular openings with a 1 mm radius, and a ME-coated sponge) field trap was hung in a tree (1.5-4 m above ground level) at the centre of the selected area.Trapped flies were collected once a month from June 2020 to December 2022, replacing new ME-coated sponges in each sampling round.Collected fruit flies were put in transparent polythene bags and then brought to the laboratory in the Department of Zoology, University of Ruhuna, for identification.These specimens were identified using taxonomic keys (Prabhakar et al. 2012;Schutze 2012;Choudhary et al. 2014;Plant Health Australia 2018;Daud et al. 2020;Leblanc et al. 2021).

Data analysis (i) Density of B. dorsalis
The density of recorded B. dorsalis was calculated using the following formula: where, D = density, l = number of specimens of B. dorsalis, L = number of all fruit fly specimens.(ii) Spatiotemporal maps The mean density of B. dorsalis in each bioclimatic zone for each study year was calculated.Then the density variations of B. dorsalis among bioclimatic zones for the years 2020, 2021, 2022, 2023, and 2024 were depicted by GIS maps using QGIS software (version 3.28.7,Firenze).Further, the trend of the density of B. dorsalis in each bio-climatic zone during the period from 2020 to 2024 was illustrated using colour intensity maps.The intensity of B. dorsalis density in the particular study year was depicted in graphs using "low, moderate, high, and very high" categories.The following four stages were followed: Data preparation; plotting the density to recognize the seasonality and data proper transformations.Model selection; computing autocorrelation function (ACF) and partial autocorrelation function (PACF), and examining and plotting.Estimation and diagnostics; estimating the pattern of autocorrelations and partial autocorrelations of moving average (MA) or autoregressive (AR) based on the ACF and PACF values, based on the significance level of autocorrelations or partial autocorrelations of AR or MA, the model was considered appropriate.Forecasting; using the selected best-fit model.Data analysis was conducted using R (version 4.0) software.Seasonality plots, ACF and PACF plots, and residual error plots for the density of B. dorsalis in each bio-climatic zone were illustrated using this software.
(iv) Checking the seasonality of B. dorsalis density Wet zone: Fig. 2 depicts a time series fluctuation plot of B. dorsalis monthly density in the wet zone from July 2020 to December 2022.The time series contains both seasonal and trend components.As a result, log-transformed data was created and used for further analysis.
Intermediate zone: Fig. 3 shows that the variance of the time series does not vary significantly, but an indistinct trend appears to begin at the end of May in 2021.Therefore, auto-correlation function (ACF) and partial auto-correlation function (PACF) plots are evaluated using log-transformed data.
Dry zone: Fig. 4 indicates that the data set has a seasonal element and a trend variation.Hence, a seasonal ARIMA model can be used to model the data.As a result, the ACF and PACF plots are first used to evaluate the model.
Arid zone: The data can be modelled using a seasonal ARIMA model since the time series plot in Fig. 5 indicates that there is a trend variation and a seasonal component.

(v) Model selection, estimation and diagnostic
Wet zone: Fig. 6 depicts the equivalent ACF and PACF.The ACF plot indicates that there is a seasonal component.There is a significant cutoff in the correlation for the PACF figure at lag 1.This means that q should be set to 1 and that the series uses the moving average (MA) (1) procedure.As a result, the seasonal ARIMA model can be used to model this series.
Intermediate zone: The plotted Fig. 7 shows the sample ACF and PACF values.Seasonality is shown in the ACF graph, and the PACF plot suggests that auto-regressive (AR) (1) might be a suitable mode for the data.As a result, the seasonal ARIMA model can be used to model this series.
Dry zone: Fig. 8 shows the derived sample ACF and PACF values.The proper MA (1) model, with q = 1, is represented by the PACF graphic.As a result, the ARIMA model can be used to represent this series.
Arid zone: The seasonality component might be taken into account based on the estimated sample ACF and PACF values that are presented in Fig. 9, and the PACF plot suggests an AR (1) model.
The best-fit seasonal ARIMA model to forecast the density in each bio-climatic zone was auto-generated by the software.The ARIMA model consists of three main components, such as auto-regression (AR), moving average (MA), and integration (I) terms.This model is useful for modelling both non-seasonal (p, d, q) and a wide range of seasonal data (P, D, Q) as well as an autoregressive term denoted as AR (p), the order of differences to make the non-stationary time series stationary as (d), and the number of moving average terms as MA (q).The fitted equation for each identified best ARIMA model was created using ARIMA notations as follows: The representation of multiplicative seasonal ARIMA as follows, where, Φ p (B s ) is the seasonal AR operator of order P; Φ p is the regular AR operator of order p; (1 − B) d X t represents the seasonal differences and d = (1 − B)d the regular differences; ϴ Q (B s ) is the seasonal moving average operator of order Q; Ѳ q (B) is the regular moving average operator of order q; e t is a white noise process.

Results
The forecasting of the density of B. dorsalis are shown for each wet, intermediate, dry, and arid zone as follows:

Best fit models and residual values
According to the preceding estimations, the auto-generated best fit models for B. dorsalis density in the bioclimatic zones are indicated in Table 2.
Further residual plots imply that these models are fair enough in forecasting the density of B. dorsalis in wet, intermediate, dry and arid zones.

Forecasting of B. dorsalis density
Wet zone: Fig. 10   Arid zone: Compared to the other three climate zones, the arid zone anticipates the greatest increase in B. dorsalis density.The density of B. dorsalis is on increasing, and in 2025, this upward tendency could reach 37% (Fig. 13).
Note: The colour shades on both sides of the trend curve in Figs. 10,11,12 and 13 indicate the level of the variation of standard deviation.

Discussion
The basic data analysis demonstrated that the seasonal ARIMA model is the best fit model for the data set since the density of B. dorsalis in all bioclimatic zones implied a seasonal component in every study year.As shown by Cooray (2008), later probability plots and residual values of each chosen model for the relevant bioclimatic zone showed the best fit for the specific data set.De Villiers et al. (2016) emphasize the need to take into account the seasonal element when calculating the density of B. dorsalis.They pointed out that studies frequently refer to population survival using occurrences and disregard seasonal change.De Villiers et al. (2016), contend that B. dorsalis exhibits a seasonal population and that identifying their advantageous seasons is essential to applying control methods to avoid significant losses to fruits, provide further support for the application of seasonal ARIMA for the current data series.
Furthermore, Dong et al. (2022) noted that because B. dorsalis distribution is strongly related to the seasonal conditions in the environment, investigations on its forecasts should concentrate on both its seasonal and yearround distribution.According to Peris (2016), the main mango fruiting season corresponds with the wet period of in each bioclimatic zone.As such, the seasonality of B. dorsalis density in each bioclimatic zone is mostly determined by the fruiting season and the wet period.Chen et al. (1995), Lv et al. (2008), Patel et al. (2013), Bana et al. (2017) and Momen et al. (2022), all support the current finding by revealing that the seasonal variation of B. dorsalis closely correlates with the harvesting periods of their main host plant.
The arid zone is expected to have the highest B. dorsalis density in 2025, indicating a high risk to the zone's future fruit sector.The wet zone had the lowest expected density of B. dorsalis in 2025, which suggests that the probability of risk is decreasing for wet zone fruits in the future.It further shows the possibility of B. dorsalis density variations over time among Sri Lanka's bioclimatic zones between 2020 and 2025.Intriguingly, in the arid zone, B. dorsalis exhibited the lowest density in 2020, and by 2025, its highest density may be predicted.Despite a notable increase in B. dorsalis density being observed in the intermediate zone for the years 2020-2022, their predicted density from 2023-2025 stays moderate.Interestingly, the population of B. dorsalis in the wet zone decrease noticeably from 2020 to 2025.These findings may be explained by significant behavioral characteristics of B. dorsalis, including its propensity for rapid invasion, broad distribution, and potential environmental climate tolerance (Loomans et al. 2019).
This foremost, but valuable, study results make it important for agriculture authorities and farmers in Sri Lanka to take necessary precautions to control the B. dorsalis populations to avoid serious economic loss in the future fruit industry in Sri Lanka.This astounding findings can be further useful for designing, planning, and implementing proper scientific control measures in managing future B. dorsalis devastating fruit damages in the bioclimatic zones of Sri Lanka.

Conclusions
Bactrocera dorsalis shows a year-round density variation with a seasonal component.B. dorsalis would have an upward density trend from 2020 to 2025.The forecast density increase of B. dorsalis is 20%, 30%, 26%, and 37%, respectively, in the wet, intermediate, dry, and arid zones.The highest B. dorsalis density from the arid zone and the lowest from the wet zone could be predicted in 2025.

Fig. 1
Fig. 1 Selected forty study sites in four bioclimatic zones of Sri Lanka (modified the basic map derived from Alahacoon and Edirisinghe 2021)

Fig. 6 Fig. 7 Fig. 8 Fig. 9
Fig. 6 Estimated sample ACF and PACF values (wet zone) shows the B. dorsalis predicting density using the aforementioned model equation for the years 2023-2025.The density of B. dorsalis is showing an upward trend from 2020 to 2025.In comparison to 2020, a 20% increase in B. dorsalis density is expected in 2025.Intermediate zone: Fig. 11 depicts a rising trend in density over time for the mean annual density of B. dorsalis from 2020 to 2025.In comparison to 2020, B. dorsalis' density could grow by 30% in 2025.Dry zone: Fig. 12 shows an increasing trend in the mean annual density of B. dorsalis from 2020 to 2025.In comparison to 2020, B. dorsalis density increases by 26% in 2025.

From
2020 to 2025, B. dorsalis density could increase in all bioclimatic zones.The intermediate zone from 2020 to 2024 has the maximum density of B. dorsalis, as seen in Fig. 14.The arid zone is expected to have the highest density between 2024 and 2025.

Fig. 14
Fig. 14 The overall variation of B. dorsalis mean density among bioclimatic zones from 2020 to 2025

Fig. 16
Fig. 16 Mean density variation of B. dorsalis from 2020 to 2025 in dry and arid zones (D: dry; A: arid)

Table 1
Description of sampling localities

Table 2
Best fit models for forecasting the B. dorsalis density in bioclimatic zones