Jing Wan,QI Guo-jun,Ma Jun,Yonglin REn,WanG Rui,Simon MCKIRDY

1 Harry Butler Institute,Murdoch University,Perth,WA 6150,Australia

2 Guangdong Provincial Key Laboratory of High Technology for Plant Protection,Plant Protection Research Institute,Guangdong Academy of Agricultural Sciences,Guangzhou 510640,P.R.China

3 Guangdong Inspection and Quarantine Technology Center,Guangzhou 510623,P.R.China

4 Institute of Plant Protection,Chinese Academy of Agricultural Sciences,Beijing 100193,P.R.China

abstract Bactrocera bryoniae and Bactrocera neohumeralis are highly destructive and major biosecurity/quarantine pests of fruit and vegetable in the tropical and subtropical regions in the South Pacific and Australia.Although these pests have not established in China,precautions must be taken due to their highly destructive nature.Thus,we predicted the potential geographic distribution of B.bryoniae and B.neohumeralis across the world and in particular China by ecological niche modeling of the Maximum Entropy (MaxEnt) model with the occurrence records of these two species.Bactrocera bryoniae and B.neohumeralis exhibit similar potential geographic distribution ranges across the world and in China,and each species was predicted to be able to distribute to over 20% of the globe.Globally,the potential geographic distribution ranges for these two fruit fly species included southern Asia,the central and the southeast coast of Africa,southern North America,northern and central South America,and Australia.While within China,most of the southern Yangtze River area was found suitable for these species.Notably,southern China was considered to have the highest risk of B.bryoniae and B.neohumeralis invasions.Our study identifies the regions at high risk for potential establishment of B.bryoniae and B.neohumeralis in the world and in particular China,and informs the development of inspection and biosecurity/quarantine measures to prevent and control their invasions.

Keywords:fruit fly,Bactrocera bryoniae, Bactrocera neohumeralis, ecological niche modeling,MaxEnt,potential geographic distribution,habitat suitability

1.Introduction

Invasive species not only pose a threat to global biodiversity but also affect the economy,natural resources,and human health (Pimentelet al.2005;Earlyet al.2016;Wan and Yang 2016).Due to the increase in global trade and transportation in the past few years,the risk and rate of invasive alien species (IAS) introductions have increased significantly (Meyerson and Mooney 2007;Westphalet al.2008;Linet al.2011;Aguin-Pombo 2012).Controlling and eradicating IAS can be challenging and costly once they are established in a new region.Predicting the potential geographical distribution of IAS before their introduction can help develop proactive measures for prevention (Hobbs and Humphries 1995;Waage and Reaser 2001;Liebholdet al.2015).Therefore,identifying locations where IAS are likely to establish is critical for the early detection,prevention,and control of IAS invasions to mitigate ecological and economic damages (Jiménez-Valverdeet al.2011;Zhuet al.2017).

Recently,ecological niche modeling (ENM) has emerged as a powerful tool for studying biological invasion and predicting the potential global and regional distribution,habitat suitability,and possible dispersal routes of IAS(Peterson 2001;Barveet al.2011;Jiménez-Valverdeet al.2011;Zhuet al.2017).Many modeling approaches,such as the Generalized Linear Model (GLM),Generalized Additive Model (GAM),Bioclimatic Prediction System (BIOCLIM),Climate Change Experiment (CLIMEX),Genetic Algorithm for Rule-Set Production (GARP),and Maximum Entropy(MaxEnt) have been employed in biodiversity,conservation biology and invasion biology studies (Bellardet al.2012;Zhu and Qiao 2016).Among these,MaxEnt is a generalpurpose machine learning method with many advantages,such as using continuous and categorical data,incorporating interactions between different variables (Phillipset al.2006),and avoiding commission errors (Pearsonet al.2007).It also outcompetes other ecological niche models with a better performance in predicting species distribution (Elithet al.2006;Hernandezet al.2006).

The genusBactroceracontains over 500 fruit fly species and many are considered major pests,leading to quantitative and qualitative losses in crop production in Australia,Oceania and tropical Asia (White and Elson-Harris 1992;Doriet al.1993;Vijaysegaran 1997).For instance,the Queensland fruit flyBactroceratryoniis one of the most destructive insect pests of fruit and vegetable crops in Australia (Raghuet al.2000;Sutherstet al.2000;Clarkeet al.2011),causing an estimated annual lose between 28.5 million and 100 million AUD (Sutherstet al.2000).Among all pest fruit fly species,B.bryoniaeandB.neohumeralisare economically significant due to wide climatic and host range(Raghu and Clarke 2001;Meats 2006;Morrowet al.2015;Royer 2015),withB.tryonicausing more damage to fruit and vegetable production (Meats 2006;Royer and Hancock 2012).Bactrocera bryoniaeandB.neohumeralisare widely distributed throughout the tropical and sub-tropical regions of Australia and South Pacific regions (Drewet al.1982;Gillespie 2003),where they either directly damage fruit and vegetable crops or cause phytosanitary restrictions of the entire infected planting area (Dominiak and Worsley 2016).

Bactrocera bryoniaeandB.neohumeralisare biosecurity/quarantine fruit fly pests that can be spread from the tropical and subtropical regions of Australia (Dominiak and Worsley 2016).They were first reported and intercepted in China in early 1998,and arrived with pepperviaa passenger from Australia (Lianget al.1998).AlthoughB.bryoniaeandB.neohumeralishave not yet colonized naturally in China,they are of major concern owing to their destructive nature and frequent interceptions at Chinese ports (Lianget al.2008;Wu 2014).As agricultural trade between Australia and China continue to increase,Australia has become the fourth largest trading partner of China for agricultural products(Yin and Xiao 2007).As a result,there is a growing risk ofB.bryoniaeandB.neohumeralisinvasions.Therefore,inspection and quarantine are required when importing fruits from Australia.

The invasion ofB.bryoniaeandB.neohumeraliswill threaten the agriculture production and devastate the fruit industry in China (Lianget al.2008;Liet al.2009).Understanding the habitat suitability forB.bryoniaeandB.neohumeralisestablishment across various geographical regions in China is crucial to prevent and control their spread (Liet al.2009).Thus,we aim to predict the potential geographical distributions of these two species based on habitat suitability.We employ MaxEnt for the prediction using global-scale occurrence records ofB.bryoniaeandB.neohumeralisaccording to the ecological niche requirements.Regions that are environmentally suitable for these two fruit fly species will be identified.The findings of this study will guide the development of the control measures for preventingB.bryoniaeandB.neohumeralisinvasions.

2.Materials and methods

2.1.Species occurrence data

The occurrence data ofB.bryoniaeandB.neohumeraliswere derived from the Global Biodiversity Information Facility(http://www.gbif.org),Center for Agriculture and Bioscience International,Invasive Species Compendium (http://www.cabi.org/isc),European and Mediterranean Plant Protection Organization,Plant Quarantine Data Retrieval system(http://www.eppo.int) and published studies.These datasets provide the most complete and detailed historical and distributional information ofB.bryoniaeandB.neohumeralisto date.

Niche model predictions tend to overfit around known occurrences,and the performance values of the model are often inflated when using spatial cluster occurrences(Jiménez-Valverdeet al.2011;Boriaet al.2014).To eliminate spatial autocorrelation and sampling bias,we created a grid of 10 km×10 km cells using ArcGIS and randomly selected one or more single point(s) from each cell.As shown in Fig.1,23 locations ofB.bryoniaeand 33 locations ofB.neohumeraliswere retained after the filtering.

2.2.Climatic data

Climatic factors are more important than topographical factors for determining the potential distribution of invasive species on a larger geographical scale (Guisanet al.2013).Non-climatic factors such as bionomics and occupancy dynamics,may function on a local scale (Hortalet al.2010).Therefore,climatic factors are widely used to predict the potential distribution ranges of invasive species on the regional and global scales.

Climatic variables used for modeling were obtained from the WorldClim Database (http://www.worldclim.org/).Nineteen climate variables recorded between 1950 and 2000 were used including annual mean temperature (Bio1),monthly mean temperature range(Bio2),isothermality (Bio3),temperature seasonality(Bio4),maximum temperature of the warmest month(Bio5),minimum temperature of the coldest month (Bio6),annual temperature range (Bio7),mean temperature of the wettest quarter (Bio8),mean temperature of the driest quarter (Bio9),mean temperature of the warmest quarter(Bio10),mean temperature of the coldest quarter (Bio11),annual precipitation (Bio12),precipitation of the wettest month (Bio13),precipitation of the driest month (Bio14),precipitation seasonality (Bio15),precipitation of the wettest quarter (Bio16),precipitation of the driest quarter(Bio17),precipitation of the warmest quarter (Bio18),and precipitation of the coldest quarter (Bio19).All climatic variables had a spatial resolution of 2.5 arcminutes (～ 5 km at the equator) (Zhuet al.2012).

2.3.Variable selection and statistical analysis

Fig.1 Global occurrence data of Bactrocera bryoniae (Black spot) and Bactrocera neohumeralis (Red spot) for the current modeling.

Multicollinearity among the predictor variables may hamper the analysis of the relationship between species distribution and environment.Therefore,some causal variables must be discarded if others can better statistically explain the variation in the response variable (Heikkinenet al.2006).To establish a high-performance model with fewer variables,cross-correlations (Pearson correlation coefficient,r) among the climatic variables were tested.

We also assessed the correlations amongB.bryoniaeandB.neohumeralisoccurrence points and each climatic variable using ArcGIS 10.4.Only one variable from each set of the highly cross-correlated variables (r>0.8) was retained for further study based on both correlation analyses (Zhuet al.2012,2017).Seven variables were included for theB.bryoniaeprediction,including Bio1,Bio2,Bio4,Bio12,Bio14,Bio15 and Bio18,and six were retained for theB.neohumeralismodeling,Bio1,Bio2,Bio7,Bio12,Bio14 and Bio15.

2.4.MaxEnt modeling

We used MaxEnt (version 3.4.1,available from http://biodiversityinformatics.amnh.org/open_source/maxent/)to predict the potential distribution ofB.bryoniaeandB.neohumeralisglobally and in China.The parameters of MaxEnt were set as follows:the available feature was auto feature including linear,quadratic,product,threshold and hinge (Phillipset al.2006;Elithet al.2011).The logistic output was used in MaxEnt,which generated a continuous map with an estimated probability of presence between 0 and 1 (Elithet al.2006;Zhuet al.2017).

The selected climatic variables that represent the main factors are considered important determinants of the distribution ofB.bryoniaeandB.neohumeralisin China and the world.The MaxEnt analysis randomly selected 75% data for model training,and the remaining 25% data were used to test the model (Zhuet al.2012;Qiet al.2015).In the MaxEnt model,the significance of the variables contributing toB.bryoniaeandB.neohumeralisdistribution is represented by the average values of the area under the curve (AUC) of ten model iterations (Elithet al.2010).This suitability analysis provides sensitivity and specificity scores,the average contribution rate of each variable to the model,and the jackknife analysis of the contribution of each variable to the model (Phillipset al.2006;Suwannatraiet al.2017).The habitat suitability forB.bryoniaeandB.neohumeraliswas reclassified into four different levels based on the predicted suitable probability:unsuitable (<0.05),low habitat suitability (0.05–0.1),moderate habitat suitability (0.1–0.2),and high habitat suitability (>0.2).

2.5.Model evaluation

The AUC of the receiver operating characteristic (ROC)generated by MaxEnt was used to estimate model performance (Wanget al.2007;Ward 2007).AUC values range from 0 to 1,AUC<0.5 suggests random prediction,an AUC value between 0.5 and 0.7 indicates poor performance,an AUC value between 0.7 and 0.9 indicates moderate performance,and AUC>0.9 indicates high performance(Swets 1988;Phillipset al.2006).

Fig.2 Receiver operating characteristic (ROC) curve verification of the prediction for Bactrocera bryoniae (A) and Bactrocera neohumeralis (B) by MaxEnt.AUC,area under the curve.

3.Results

3.1.Model performance for potential distribution and environmental variables

Our models showed high predictive accuracy with AUC values above 0.9,indicating that they perform better than random models in predicting the habitat suitability for the two fly species.Specifically,we detect a mean training AUC value of 0.98 and 0.99 (Fig.2-A),and an average test AUC value of 0.99 and 0.91 (Fig.2-B) forB.bryoniaeandB.neohumeralis,respectively.

According to the relative contribution of each selected climatic variable in predicting the geographic distribution of the two fly species (Fig.3-A and B),precipitation of the warmest quarter (Bio18;68.3%) was the most important variable in explaining the potential distribution ofB.bryoniae(Fig.3-A),while Bio4,Bio1,and Bio14 contributed 27.1% to the prediction,these variables contributed a total of 95.4% to the performance of the model.Annual precipitation (Bio12;46.2%) was the most important environmental variable in explaining the potential distribution ofB.neohumeralis(Fig.3-B),while Bio1,Bio14,and Bio7 contributed 48.7%to model performance,these four factors explained 94.9%of the prediction.

3.2.The global distribution of B.bryoniae and B.neohumeralis

Fig.3 The contribution of each environmental variable in predicting the habitat suitability of Bactrocera bryoniae (A) and Bactrocera neohumeralis (B) in the models.Bio18,precipitation of the warmest quarter;Bio4,temperature seasonality;Bio1,annual mean temperature;Bio14,precipitation of the driest month,Bio15,precipitation seasonality;Bio12,annual precipitation;Bio2,monthly mean temperature range;Bio7,annual temperature range.

The global projection results showed that Asia,Africa,South America,North America,and Oceania had a wide range of habitable regions forB.bryoniae(Fig.4),occupying 23.04%of the global land surface area,among which 12.29,5.24 and 5.51% had high,moderate and low suitability,respectively.The highly suitable habitats forB.bryoniaewere mainly distributed in southern Asia,the equatorial region and the southeastern coastal areas of Africa,southern North America,northern and central South America,and eastern coastal areas of Australia.The Mediterranean coastal areas,northern and western Australia,southern Africa,central North America,and parts of South America exhibited moderate to low suitability forB.bryoniae.

Compared toB.bryoniae,B.neohumeralishad a relatively wider distribution range (Fig.5).The total suitable area ofB.neohumeralisaccounted for 25.66% of the global land surface area,of which 15.27,5.80 and 4.59% were high,moderate and low suitability areas,respectively.Despite an overall wider distribution range,the highly suitable habitats forB.neohumeraliswere similar to those ofB.bryoniae,we detected no significant difference in the moderate and low suitable habitats betweenB.bryoniaeandB.neohumeralis.

3.3.Potential distribution of B.bryoniae and B.neohumeralis in China

The projection results showed that south of the Yangtze River region had a wide range of habitable areas forB.bryoniae(Fig.6).The total suitable area ofB.bryoniaeaccounted for 22.06% of the total land area of China.Of the total suitable area ofB.bryoniae,11.81,5.48 and 4.77% had high,moderate and low suitability,respectively.The highly suitable areas forB.bryoniaeare predominantly distributed in Hainan,Guangdong,Guangxi,Fujian and Taiwan,most parts of Yunnan,southern Jiangxi,Guizhou and Zhejiang,southeastern Sichuan,and some scattered areas in Tibet Autonomous Region.By contrast,the northern Yangtze River region had the smallest habitat forB.bryoniae(Fig.6).

Fig.4 Global habitat suitability map for Bactrocera bryoniae.

Fig.5 Global habitat suitability map for Bactrocera neohumeralis.

Compared toB.bryoniae,the suitable habitat range forB.neohumeraliswas slightly larger (Fig.7).The total suitable area ofB.neohumeralisaccounted for 22.29% of the total land area of China,of which 14.87,5.44 and 1.98%had high,moderate and low suitability,respectively.The highly suitable areas forB.bryoniaewere mainly distributed in Hainan,Guangdong,Guangxi,Fujian,and Taiwan,most parts of Jiangxi,Hunan,Guizhou,Yunnan,eastern Sichuan,most parts of Chongqing,Southern Hubei,eastern Zhejiang,and some scattered areas in Jiangsu,Shanghai and Tibet(Fig.7).Similar toB.bryoniae,northern Yangtze River had the least suitable area forB.neohumeralis(Fig.7).

Fig.6 Potential habitat suitability map for Bactrocera bryoniae in China.

Fig.7 Potential habitat suitability map of Bactrocera neohumeralis in China.

4.Discussion

MaxEnt utilizes continuous and categorical data,incorporates interactions between different variables (Phillipset al.2006),avoids commission errors (Pearsonet al.2007),and outperforms other ecological niche models in predicting species distribution (Elithet al.2006;Hernandezet al.2006).We used MaxEnt to predict the potential geographic distribution ofB.bryoniaeandB.neohumeralisdue to the limited availability of sample localities and the lack of biological characteristic information such as temperature tolerance range.In this study,the MaxEnt model performed well in predicting suitable habitats forB.bryoniaeandB.neohumeraliswith AUC values above 0.9 (Fig.2),which is referred to as an excellent model-fit category (Phillips and Dudík 2008).

Given their destructive impacts on local commercial fruit and vegetable crops in Australia (Morrowet al.2015;Royer 2015),the risk ofB.bryoniaeandB.neohumeralisinvasions in the rest of the world must be considered (Zaluckiet al.1984;Meats 2006;De Meyeret al.2008).Upon invasion,an invasive species must overcome a series of challenges such as geographic barrier,climate,and natural enemies(Richardson and van Wilgen 2004;De Meyeret al.2010).According to our models,regions with high suitability are more prone to have invasion outbreaks of these fly species.Specifically,southern Asia,the central and southeastern coast of Africa,southern North America,northern and central South America,and the eastern coast of Australia are predicted as highly suitable regions for bothB.bryoniaeandB.neohumeralis(Figs.4 and 5).It is worth noting that the native habitats ofB.bryoniaeandB.neohumeralisin Australia were among the predicted invasion areas,reflecting the detection power of our models.Despite the high detection power,this study cannot eliminate the possibility ofB.bryoniaeandB.neohumeralisinvasions in regions shown negatively in our prediction.

The results of MaxEnt model revealed similar distribution patterns ofB.bryoniaeandB.neohumeralisin China.However,B.neohumeralishas a wider suitable habitat range comparing withB.bryoniae(Figs.6 and 7).In China,most areas of the southern Yangtze River are suitable habitat forB.bryoniaeandB.neohumeralis.This finding is similar to the potential geographical distribution of the Queensland fruit fly (B.tryoni) (Raoet al.2009).Meats (2006) reported a widespread distribution ofB.bryoniaeandB.neohumeralisin tropical rainforest and subtropical monsoon warm climate in the eastern coast of Australia.Thus,increased habitat suitability of these two fly species across the southern Yangtze River area is likely owing to its similar climatic condition to that of the Queensland fruit fly (Yan 2006;Raoet al.2009).Moreover,southern China has both tropical and subtropical climates and produces a variety of tropical fruits (Qiet al.2015),making it the most vulnerable target ofB.bryoniaeandB.neohumeralisinvasions.

The accuracy of the model depends on the quality and quantity of the occurrence records and the selection of environmental variables (Jiménez-Valverdeet al.2011;Qiaoet al.2013;Zhu and Qiao 2016).Recent studies suggested that temperature and precipitation have a significant impact on the physiology of fruit fly species (Guisanet al.2013).But Meats (2006) argued that winter field temperatures do not limit the geographical range ofB.neohumeralis.In Queensland,Australia,B.neohumeralishas been reported in areas where the average annual rainfall exceeds 760 mm(Drew 1982).In this study,the MaxEnt model takes into account of temperature,precipitation,and other variables that potentially influence the distribution ofB.bryoniaeandB.neohumeralis.Precipitation of the warmest quarter and annual precipitation were the most important among all environmental variables,contributing 68.3 and 46.2% to the prediction ofB.bryoniaeandB.neohumeralisdistribution,respectively (Fig.3).Consistent with this notion,subtropical regions in South and Southwest China,where rainfall is sufficient,were predicted to have an increased fruit fly occurrence.According to these results,precipitation does not seem to be a limiting factor ofB.bryoniaeandB.neohumeralisoccurrence in South and Southwest China.

Fruit fly pests pose a great threat to fruit and vegetable production throughout the world,causing both quantitative and qualitative losses (White and Elson-Harris 1992;Vijaysegaran 1997;Leblancet al.2012).The most common hosts ofB.bryoniaeinclude pepper,mango,passionflower,and banana (Lianget al.1998),while similar toB.tryoni,hosts ofB.neohumeralisinclude a wide range of plant species and is highly destructive (Drew 1982;Allwoodet al.1999).Consistently,our results show thatB.neohumeralishas a higher invasion risk thanB.bryoniaeby showing a broader host range and a higher climate tolerance.To preventB.neohumeralisandB.bryoniaeinvasion,strict inspection and biosecurity/quarantine measures must be enforced at the borders,especially when importing fruits from Australia,New Zealand and the South Pacific regions.

5.Conclusion

The results of this study reveal a similar distribution pattern ofB.bryoniaeandB.neohumeralisin China and the world.However,in China,particularly the southern Yangtze River area was predicted to be the most suitable forB.bryoniaeandB.neohumeralissurvival.Therefore,South China was considered to have the highest risk of invasion.Compared with the potential geographic distribution,the host adaptability and climate tolerance ofB.bryoniae,B.neohumeralismay have high a risk of invasion into these predicted regions.

Acknowledgements

This research was supported by the National Key R&D Program of China (2017YFC1200600 and 2016YFC1202104),the Innovation Team of Modern Agricultural Industry Generic Key Technology R&D of Guangdong Province,China(2019KJ134),and the Open Fund of the Guangxi Key Laboratory of Biology for Crop Diseases and Insect Pests,China (2016-KF-3).A student scholarship was provided by the Harry Butler Institute,Murdoch University,Australia.