Mapping the Potential Global Range of the Brown Marmorated Stink Bug , Halyomorpha halys , with Particular Reference to New Zealand

Originating from Asia, the brown marmorated stink bug (BMSB) is a significant pest of horticultural/agricultural crops, grapes, woody ornamental and herbaceous plants, and is also a nuisance to people, due to its overwintering behavior in human habitation. The global range of this pest is steadily increasing and previous predictions of environmental suitability have shown New Zealand to be highly suitable. Due to the economic value of horticultural and agricultural industries to the New Zealand economy, it is vital to understand the range of potential risk within the country. Global and New Zealand potential suitability for BMSB was modeled using three algorithms and the resulting predictions ensembled to predict the potential range under current climatic conditions and under trajectories of future low (Representative Concentration Pathways, RCP, 2.6) and high (RCP 8.5) greenhouse gas emissions for both 2050 and 2070. Under current conditions, models showed a high global suitability within latitudes 25◦–50◦ N, southern South America, southeast and southwest regions of Australia and large areas of New Zealand. Modeling the effect of climate change on BMSB range in New Zealand resulted in a southerly range shift over time, particularly with high emissions trajectory. Currently, BMSB is not established in New Zealand and it is vital that this remains the case.


Introduction
New Zealand, being an island nation heavily reliant on primary industry productivity, is particularly vulnerable to the impact of invasive species.Currently, New Zealand's Ministry of Primary Industries (MPI) is on high alert to prevent the entry of Halyomorpha halys or brown marmorated stink bug (BMSB) from entering the country [1].This arthropod, of approximately 17 mm in length [2], is a native of Japan, Korea, China and Taiwan [3] and is known as a horticultural pest and nuisance to humans in an increasing number of locations around the world.First introduced to Allentown, Pennsylvania in the 1990s [3], it is now found in 43 states in the USA and four Canadian provinces [4].The origin of the BMSB in the US has been traced to the Beijing area of China [5].Although older records are available for presence of this arthropod in Lichtenstein in 2004 [6], the BMSB was first officially reported in 2008 from collections in Europe in the Zürich region of Switzerland in 2007 [7].Subsequently, it has spread within Switzerland, was recorded in Germany in 2012 [8], Strasbourg in France in 2012/2013 [9], northern Italy in 2013 [10] and Hungary in 2014 [11].Detection of BMSB was recorded in the Dunedin area, New Zealand, in 2010, [12] while reported interceptions by MPI have shown steady increases in the 2014 to 2016 period [13].
As good fliers and generalist feeders on a wide range of plants, such as vegetables, fruit and ornamental plants [14], this species has the potential to easily disperse and establish in novel environments.Despite this high association with plant material, it is suggested that entry to New Zealand is less likely on fruit and nursery stock, due to the high health standards applied to imported goods [14].Entry to New Zealand is, however, more likely due to the species overwintering behavior where it infiltrates personal or commercial items such as luggage, freight and loaded containers [14].
Known as a destructive pest of various fruits, particularly apple trees, and soybean crops in its native range [15], the brown marmorated stink bug has proven to be an increasingly significant pest of a range of orchard crops, agricultural crops such as soybeans and sweet corn, grapes and ornamental plants [16,17].Damage is varied but may include scarring, depressions and internal tissue damage, premature flower drop, reduction in field corn kernel size and color change, "stay green" syndrome in soybeans and early fruiting body drop, resulting in low yield in vegetable and fruit crops, and potential tainting of grape juice [16].This damage can result in significant economic loss, as reported in the 2010 apple crops in the mid-Atlantic regions of the US, which was estimated at approximately $37 million [18], and product losses exceeding 90% in some stone fruit crops [19].The range of plant host species, including woody and herbaceous plants, appears to be steadily increasing [20], with stippling and scabbing leaf damage from foliage feeding by the arthropod [3] to profuse sap flow from injury sites due to bark feeding [21].In addition to fruit and crop damage, BMSB is also known to be a nuisance to people, due to its overwintering behavior [22].The bug tends to swarm together during winter in human habitation, and produces of an unpleasant odor if crushed [22][23][24].
Ecological Niche Modeling (ENM) has been used to determine potential range distribution of invasive species in new environments [25][26][27].These predictions provide valuable information for the protection of novel environments and the prevention and/or reduction of impacts resulting from invasive species.The potential global invasive range of the BMSB has previously been mapped by Zhu et al. [28].The authors predicted that northern Europe, northeastern North America, Angola, Uruguay, southern Australia and the North Island of New Zealand showed high climate suitability.
Due to the potential severity of the impact of this species both in terms of plant susceptibility, economic loss, a nuisance to human society, and the increase in New Zealand border detection of BMSB, it is important to determine the potential suitable environmental range in New Zealand.This situation is of greater concern under the influence of climate change, whereby the environmental suitable is expected to increase across the country [29][30][31][32].This study aims to determine the potential global and New Zealand distribution of BMSB under current climatic conditions and future environmental suitability under the influence of two climate change scenarios; low (Representative Concentration Pathways, RCP, 2.6) and high (RCP 8.5) greenhouse gas emissions for both 2050 and 2070.This will highlight the extent of the potential problem, indicate the risk of invasion, and support the management for prevention of entry, as well as its future regional and national management if this destructive species become established in New Zealand.

Materials and Methods
Occurrence data for this species in both its native and invasive range was collected using Global Biodiversity Information Facility (http://www.gbif.org/)and the US Department of Agriculture Animal and Plant Health Inspection Services [33].Occurrence points were checked for redundancy and spatially rarified using the SDMtoolboxv1.1 [34] in ArcMap 10.2.A total of 573 records were downloaded, and after manual checking and SDMToolbox correction for spatial autocorrelation and redundancy, 125 were used for modeling [34].The deletion of points determined to contribute to spatial autocorrelation was necessary to increase confidence that more robust prediction models result with less inflated performance measures are produced from spatially independent data [35][36][37].
The environmental variables used were derived from the WorldClim Database for current conditions and future climatic scenarios (http://worldclim.org/).The world model was at a resolution of 5 min and the New Zealand projection at 30 s.The future climatic scenarios consisted of the data used in the fifth Intergovernmental Panel on Climate Change (IPCC) assessment report that was based on the representative concentration pathways (RCP 2.6 representing lower emissions and RCP 8.5 representing higher emissions) [38].Models were run for both emission trajectories for the years 2050 and 2070 to predict environmental suitability.Environmental variables at the occurrence points, determined using the ArcMap tool Extract Multi-values from rasters, were used in the models (Table 1).These consisted of 19 parameters, which were checked for cross correlation using the SDMtoolbox in ArcMap [34].Only variables that were not cross-correlated at Pearson coefficient of less than or equal to 0.7 were used in the model, and with a variance inflation factor (VIF) of less than 10 [39].Studies have shown that using simpler models with a smaller number of variables improves projection of species distribution models [25,40,41].Of the 19 variables, seven (BIO1, BIO2, BIO8, BIO9, BIO12, BIO15, and BIO19) were selected in the final model (Table 1).To check for co-linearity, the VIF of the environmental variables was checked with the r-package "usdm" [42].VIF values for the non-correlated variables ranged from 3.11 to 7.06 or below commonly used thresholds [42].Species distribution modeling was run using the R package "dismo" using three high performing algorithms; Random Forest (RF), Support Vector Machine (SVM) and Maxent, using the default values for each algorithm [43,44].These have been found to be high performance algorithms with occurrence only data [45].A worldwide distribution model was generated for each of the algorithms and projected to New Zealand for current and future conditions [32,43,44].Each model was evaluated according to Area Under the Curve (AUC), Kappa and True Skills Statistic (TSS).The AUC, while some problems are reported with its use, is a threshold independent widely employed performance metric that shows the accuracy of the models in correctly predicting presence and absence (values of >0.7 showing good performance) [46,47].The "Kappa" statistic describes agreement between the observations and predictions, and incorporates the effect of chance on the accuracy of the prediction (values nearing 1.0 showing good performance) [48,49].The TSS is a commonly used measure that overcomes some of the weaknesses of the Kappa statistics (values near 1.0 indicate good performance) while addressing Kappa's sensitivity to prevalence [50].
Prediction maps of environmental suitability for each of the models were ensembled, each being weighted using the combined ranked values of AUC, Kappa and TSS [51,52].Ensembled approaches have the advantage of addressing the individual weaknesses of each modeling approach or deriving results similar to a consensus [51,52].The ensemble modeling approach has been used to predict the suitability of areas for the Argentine ant in Spain [53], invasive aquatic fishes in United States [54], pine forests in China [55], and forecasting distributions of the Queensland fruit fly in climate change scenarios [32].The R package "dismo" [56] was used with the built-in Maxent, SVM and RF algorithms to generate prediction maps and derive the final ensembled maps for each climatic scenario.The resulting ensembled predictions, in the form of rasters, were processed for presentation in ArcMap 10.2.The predictive maps were reclassified into low, medium and high categories using equal distribution categories in order to clearly present the differences between current and future suitability and measure the area that changed between periods [43].Commonly used thresholds to determine presence absence were not utilized because of the reported insensitivity of some thresholds to the lack of absence data as well as the reported variability between different thresholds [57,58].The centroid of the high category area for each projection was also determined and plotted to show the direction and magnitude of movement or shift of the highly suitable areas [32,43,44].Territorial boundaries were overlaid on the suitability model to provide information for management purposes.

Occurrence Data
Global spatially rarified occurrence points used in the model consisted of 125 sites, which were mainly located in the Northern Hemisphere (Figure 1).
Climate 2017, 5, 75 4 of 14 scenarios [32].The R package "dismo" [56] was used with the built-in Maxent, SVM and RF algorithms to generate prediction maps and derive the final ensembled maps for each climatic scenario.The resulting ensembled predictions, in the form of rasters, were processed for presentation in ArcMap 10.2.The predictive maps were reclassified into low, medium and high categories using equal distribution categories in order to clearly present the differences between current and future suitability and measure the area that changed between periods [43].Commonly used thresholds to determine presence absence were not utilized because of the reported insensitivity of some thresholds to the lack of absence data as well as the reported variability between different thresholds [57,58].
The centroid of the high category area for each projection was also determined and plotted to show the direction and magnitude of movement or shift of the highly suitable areas [32,43,44].Territorial boundaries were overlaid on the suitability model to provide information for management purposes.

Occurrence Data
Global spatially rarified occurrence points used in the model consisted of 125 sites, which were mainly located in the Northern Hemisphere (Figure 1).

Model Evaluation
Model evaluation shows that RF was the highest performer in terms of AUC (Table 2) followed closely by Maxent and SVM.Kappa and TSS is highest for RF followed by SVM.The minor differences shown with AUC reflect similar performance between the algorithms.

Model Evaluation
Model evaluation shows that RF was the highest performer in terms of AUC (Table 2) followed closely by Maxent and SVM.Kappa and TSS is highest for RF followed by SVM.The minor differences shown with AUC reflect similar performance between the algorithms.

Spatial Distribution Using Random Forest (RF) Model
The current spatial distribution model using RF (Figure 2) shows high environmental suitability in the eastern states and West Coast USA, coastal regions of southern South America, Western Europe, Japan and eastern China, the southeast and southwest coasts of Australia and considerable parts of New Zealand.The New Zealand distribution under current environmental conditions shows moderate to high environmental suitability for BMSB across the majority of New Zealand, particularly North Island, Nelson/Marlborough, coastal Canterbury and Southland regions (Figure 2).Due to the widespread environmental suitability under current conditions, there appears to be little change under predictions for low emissions trajectory for 2050.However, there is an indication of reduced environmental suitability at the top of North Island and an apparent reduced intensity of suitability over the whole of North Island associated with low emissions trajectories for 2070 and high emissions trajectories for both 2050 and 2070 (Figure 2).

Spatial Distribution Using Random Forest (RF) Model
The current spatial distribution model using RF (Figure 2) shows high environmental suitability in the eastern states and West Coast USA, coastal regions of southern South America, Western Europe, Japan and eastern China, the southeast and southwest coasts of Australia and considerable parts of New Zealand.The New Zealand distribution under current environmental conditions shows moderate to high environmental suitability for BMSB across the majority of New Zealand, particularly North Island, Nelson/Marlborough, coastal Canterbury and Southland regions (Figure 2).Due to the widespread environmental suitability under current conditions, there appears to be little change under predictions for low emissions trajectory for 2050.However, there is an indication of reduced environmental suitability at the top of North Island and an apparent reduced intensity of suitability over the whole of North Island associated with low emissions trajectories for 2070 and high emissions trajectories for both 2050 and 2070 (Figure 2).

Spatial Distribution Using Support Vector Machine (SVM) Model
The spatial distribution model using SVM for the current global distribution (Figure 3) shows similar areas of high suitability as those shown by the RF model (Figure 2), although there appears to be the presence of a diffuse pattern of moderate suitability across much of the globe with this SVM model.In New Zealand, under current conditions, high environmental suitability for BMSB is shown over a considerable proportion of North Island and a band running down the center of South Island.Projections for low and high emissions appear to show a reduction in suitability of the top of North Island and an intensification and increase in suitability in South Island over time (Figure 3).

Spatial Distribution Using Support Vector Machine (SVM) Model
The spatial distribution model using SVM for the current global distribution (Figure 3) shows similar areas of high suitability as those shown by the RF model (Figure 2), although there appears to be the presence of a diffuse pattern of moderate suitability across much of the globe with this SVM model.In New Zealand, under current conditions, high environmental suitability for BMSB is shown over a considerable proportion of North Island and a band running down the center of South Island.Projections for low and high emissions appear to show a reduction in suitability of the top of North Island and an intensification and increase in suitability in South Island over time (Figure 3).

Spatial Distribution Using Maxent Model
The spatial distribution model using Maxent for current global distribution (Figure 4) shows similar regional distribution of high suitability but a reduced global range of suitability compared with the RF and SVM models (Figures 2 and 3).Similarly, with New Zealand, this model under current conditions appears to show marked reduction on environmental suitability with only areas in South Island showing high suitability compared with the other two models (Figure 4).Addition of climate change projections appear to reduce environmental suitability in North Island and increase suitability in South Island, particularly Southland, with increasing emissions projection and time (Figure 4).

Spatial Distribution Using Maxent Model
The spatial distribution model using Maxent for current global distribution (Figure 4) shows similar regional distribution of high suitability but a reduced global range of suitability compared with the RF and SVM models (Figures 2 and 3).Similarly, with New Zealand, this model under current conditions appears to show marked reduction on environmental suitability with only areas in South Island showing high suitability compared with the other two models (Figure 4).Addition of climate change projections appear to reduce environmental suitability in North Island and increase suitability in South Island, particularly Southland, with increasing emissions projection and time (Figure 4).

Spatial Distribution Using Maxent Model
The spatial distribution model using Maxent for current global distribution (Figure 4) shows similar regional distribution of high suitability but a reduced global range of suitability compared with the RF and SVM models (Figures 2 and 3).Similarly, with New Zealand, this model under current conditions appears to show marked reduction on environmental suitability with only areas in South Island showing high suitability compared with the other two models (Figure 4).Addition of climate change projections appear to reduce environmental suitability in North Island and increase suitability in South Island, particularly Southland, with increasing emissions projection and time (Figure 4).

Ensemble Prediction
The results of models were combined to produce an ensemble model (Figure 5), which was weighted by the AUC, Kappa and TSS (Table 2).In terms of current global suitability, predominance of high suitability is shown across the globe in a latitude range of approximately 25-50 • N, east and west coasts of both southern South America and Australia, and large parts of both the North and South Islands of New Zealand.Under current climate conditions, most of the North Island appears to be highly suitable for BMSB, while the Nelson/Marlborough, Canterbury and southwest regions of the South Island show high suitability.It appears from the combined models that the suitability in North Island reduces from high to predominately moderate, while suitability increases in South Island under low and high emissions trajectory for both 2050 and 2070 (Figure 5).It appears that the impacts of climate change over time will result in an overall shift of suitable environment in a southerly direction across New Zealand.

Ensemble Prediction
The results of models were combined to produce an ensemble model (Figure 5), which was weighted by the AUC, Kappa and TSS (Table 2).In terms of current global suitability, predominance of high suitability is shown across the globe in a latitude range of approximately 25-50° N, east and west coasts of both southern South America and Australia, and large parts of both the North and South Islands of New Zealand.Under current climate conditions, most of the North Island appears to be highly suitable for BMSB, while the Nelson/Marlborough, Canterbury and southwest regions of the South Island show high suitability.It appears from the combined models that the suitability in North Island reduces from high to predominately moderate, while suitability increases in South Island under low and high emissions trajectory for both 2050 and 2070 (Figure 5).It appears that the impacts of climate change over time will result in an overall shift of suitable environment in a southerly direction across New Zealand.

Centroid Shift of Environmental Suitability for BMSB in New Zealand as Determined Using the Ensemble Model
When the high suitability category centroid was determined for all projections, there is an distinct movement to the south for all trajectories for 2050 and 2070 (Figure 6).The higher scenario RCP 8.5 showed greater magnitude of movement for the high suitability category compared to RCP 2.6.The rate of movement was estimated at 1.03 km/year for RCP 2.6 and 5.9 km/year for RCP 8.5 when measured for 2050 and 2070.

Centroid Shift of Environmental Suitability for BMSB in New Zealand as Determined Using the Ensemble Model
When the high suitability category centroid was determined for all projections, there is an distinct movement to the south for all trajectories for 2050 and 2070 (Figure 6).The higher scenario RCP 8.5 showed greater magnitude of movement for the high suitability category compared to RCP 2.6.The rate of movement was estimated at 1.03 km/year for RCP 2.6 and 5.9 km/year for RCP 8.5 when measured for 2050 and 2070.

Ensembled Model Prediction Map Combined with New Zealand Official Territories
When the ensembled predictive maps were combined with the New Zealand official territories under current conditions, they show high suitability across most of North Island territories except the upper Northland region and the Nelson/Marlbourgh, Canterbury and Southlanbd regions of South Island (Figure 7, left).Under conditions of high emissions output projected for 2070, suitable environmental range shows a southward shift in both the North and South Islands (Figure 7, right).The suitability in the Northland and Auckland regions under these conditions appears to reduce markedly.

Ensembled Model Prediction Map Combined with New Zealand Official Territories
When the ensembled predictive maps were combined with the New Zealand official territories under current conditions, they show high suitability across most of North Island territories except the upper Northland region and the Nelson/Marlbourgh, Canterbury and Southlanbd regions of South Island (Figure 7, left).Under conditions of high emissions output projected for 2070, suitable environmental range shows a southward shift in both the North and South Islands (Figure 7, right).The suitability in the Northland and Auckland regions under these conditions appears to reduce markedly.

Ensembled Model Prediction Map Combined with New Zealand Official Territories
When the ensembled predictive maps were combined with the New Zealand official territories under current conditions, they show high suitability across most of North Island territories except the upper Northland region and the Nelson/Marlbourgh, Canterbury and Southlanbd regions of South Island (Figure 7, left).Under conditions of high emissions output projected for 2070, suitable environmental range shows a southward shift in both the North and South Islands (Figure 7, right).The suitability in the Northland and Auckland regions under these conditions appears to reduce markedly.

Discussion
Modeling techniques have frequently been used to determine the spread of invasive insects under current environmental conditions [59,60] and under the influence of climate change [32,61,62].The current study has shown that the use of modeling techniques, particularly combined methods, allows the prediction of potential current and future environmental suitability for an invasive insect, BMSB, in a novel environment, i.e., New Zealand.The process of checking of occurrence data, cross correlation, co-linearity and selection of environmental variables, performance evaluation and combining of modeling methods contributes to confidence in the model outputs within this study, particularly the ensembled model.
The results of this study support the findings of Zhu et al. [28], who modeled the potential global distribution of BMSB.These authors predicted high global environmental suitability between the latitudes of 30-50 • N, which is a close match to the predictions from the ensembled model in this study (Figure 5), indicating high global suitability within latitudes 25-50 • N.Both studies predict areas of suitability in southern latitudes, including southern South America, southeast and southwest regions of Australia, and large areas of New Zealand.The potential suitability, and hence, vulnerability of New Zealand to this arthropod as indicated by Zhu et al. [28], is confirmed by the current study in the prediction of current suitability using RF, SVM, Maxent and ensembled models (Figures 2-5).
Under current climatic conditions, the models show a moderate to high general suitability for BMSB across most of the North Island and considerable areas of the South Island of New Zealand (Figures 2-5).Horticultural, orchard and viticulture production is predominant in the Auckland, Bay of Plenty, and Hawke's Bay regions in the North Island [63], all of which show moderate to high suitability for BMSB.The high producing areas of the Nelson/Marlborough and Canterbury regions of the South Island predominantly show moderate to high suitability under current climate conditions.However, future climate change predictions result in a general southerly movement of suitability for BMSB across the country (Figures 2-7), increasing the likelihood of establishment of this pest and significant impact in the future.The magnitude of shift is dependent on the emission trajectory with a higher magnitude for the higher emission trajectory.This shift in range of suitable environment is not unexpected, as species range shift has been shown to reflect predictions of climate change [29,30,32,64].
The potential negative impact of this arthropod to New Zealand's economic export earnings from fruits, vegetables, wine and seed is significant.The current export value of horticultural products (kiwifruit, wine, apples and pears, fresh and processed vegetables, other horticultural products and grains and seeds) totaled NZ$5187 million in the year ending 30 June 2016 [65].Even with conservative production losses, the severe potential economic and social impact of this organism to New Zealand is obvious.The BMSB is known to utilize a range of woody ornamental and herbaceous species [17,20] and, as such, has the potential for significant impact in the garden trade industries.
Having focused on the economic impacts of this pest, the social impacts must also be addressed.Not only will crop damage and increased pest control cost reduce profitability for producers, this arthropod is also known as a nuisance to humans, due to its overwintering behavior in cool dry refuges in human habitation, such as wall cavities, attics and garages, and its production of a characteristic unpleasant odor when crushed [17,22].Establishment of the BMSB in New Zealand will result in both economic and social impacts.
To successfully establish and spread in a novel environment, invasive species must overcome several abiotic and biotic barriers [66,67].In particular, climate/environmental matching has been an important predictor of invasive insect range expansion [31,32].Examination of the climate/environmental suitability, recorded invasion success, biological characteristics, and the impacts demonstrated by this species indicates that there is a high likelihood of establishment of the BMSB in New Zealand.Entry of BMSB to New Zealand, as predicted by Landcare Research [14], is most likely via unintentional transport in personal or commercial baggage, freight and loaded containers.Air and seaports, particularly with the current increase in visitor and immigrant numbers [68] and trade statistics [69], are, therefore, the most likely point of entry to New Zealand.Due to the characteristics of this species, such as being multivoltine (4-6 generations in southern China, 1-2 in the USA) [16], good fliers, generalist feeders [14], wide temperature tolerance (Figure 1) and the climatic suitability shown in this study, the entry at air-and seaports and the establishment of this species in New Zealand is highly likely.Furthermore, its predominant association with urban environments during establishment in new areas and likely spread by road and rail [70] support the probable ease of invasion of this species in New Zealand.As such, it can be confidently predicted that the BMSB is a significant risk to New Zealand horticultural production and as a nuisance to humans.
MPI reported that BMSB interception events in the summers of 2014/2015 and 2015/2016 were similar, but increased in winter in the 2015/2016 period [13].If these reported interception rates are calculated on an annual calendar year basis (January-December), the rates were 17 in 2014 and 65 in 2015, an increase of approximately 280%.This increase in detection may have been due to increased awareness of the threat of this organism and, therefore, increased detection effort.However, it has been reported that border interception in the summer of 2016/2017 was more than 500 incidents [71], which is significantly above detection rates in previous periods.In addition, two post-border identifications demonstrate that the arthropod is escaping border control measures [71].The resulting prediction clearly indicates that the risk of establishment of BMSB in New Zealand is increasingly alarmingly, with subsequent significant concern for the economic impacts on New Zealand agricultural and horticultural industries and potential nation-wide social impacts.
It is known that the likelihood of establishment and rate of spread of an invasive species relates to propagule pressure (propagule number and propagule size) [67,[72][73][74].Therefore, the increase in interceptions reported by MPI suggests an increased likelihood of border slippage and establishment.MPI also reported a change in origin of the interceptions with a reduction in those originating from the USA (−53%) and an increase in interceptions from Europe (+900%) and Asia (40%) in the 2014/2015 and 2015/2016 periods [13].An increase in propagule pressure, in both numbers and range of origins is likely to increase the success of establishment [72] and, therefore, increases the risk of invasion to New Zealand.
The transfer of BMSB to New Zealand native plant species is unknown.However, with evidence of an ever-increasing host species plant list [20] there is potential significant concern for the possibility that BMSB may transfer to New Zealand native woody and herbaceous plant species with unknown future impacts.If this is the case, particularly under future climate predictions, native New Zealand species may be at risk.
This modeling exercise has provided information on the environmental suitability of New Zealand under current and future projected environmental conditions.As such, it provides both national and regional authorities an indication as to the severity and impact that could be expected in their respective areas (Figure 7).This will be an important factor in the consideration of spread and management should this arthropod become established in New Zealand.This modeling exercise was focused on the effects of climate change in the New Zealand situation.However, it is realized that other important environmental variables, such as vegetation, needs to be addressed in future studies.
It is evident that New Zealand is highly susceptible to the establishment of the BMSB and, with current knowledge of its invasive impact, is extremely likely to cause significant damage to a range of productions crops, ornamental woody and herbaceous plants and potentially native species, as well as resulting in significant social impact.Due to increasing interception rates and probable border leakage, MPI is implementing significant pre-border and border strategies to minimize the risk of invasion by this species [75].However, it may be inevitable that this arthropod pest becomes established in New Zealand.The battle continues.

Figure 2 .
Figure 2. Global and New Zealand predictions of brown marmorated stink bug habitat suitability under current conditions and emissions trajectories of Representative Concentration Pathways (RCP) 2.6 and RCP 8.5 for the years 2050 and 2070 using Random Forest algorithm.

Figure 2 .
Figure 2. Global and New Zealand predictions of brown marmorated stink bug habitat suitability under current conditions and emissions trajectories of Representative Concentration Pathways (RCP) 2.6 and RCP 8.5 for the years 2050 and 2070 using Random Forest algorithm.

Figure 3 .
Figure 3. Global and New Zealand predictions of brown marmorated stink bug habitat suitability under current conditions and emissions trajectories of Representative Concentration Pathways (RCP) 2.6 and RCP 8.5 for the years 2050 and 2070 using Support Vector Machine (SVM) algorithm.

Figure 4 .
Figure 4. Global and New Zealand predictions of brown marmorated stink bug habitat suitability under current conditions and emissions trajectories of Representative Concentration Pathways (RCP) 2.6 and RCP 8.5 for the years 2050 and 2070 using Maxent algorithm.

Figure 3 .
Figure 3. Global and New Zealand predictions of brown marmorated stink bug habitat suitability under current conditions and emissions trajectories of Representative Concentration Pathways (RCP) 2.6 and RCP 8.5 for the years 2050 and 2070 using Support Vector Machine (SVM) algorithm.

Figure 3 .
Figure 3. Global and New Zealand predictions of brown marmorated stink bug habitat suitability under current conditions and emissions trajectories of Representative Concentration Pathways (RCP) 2.6 and RCP 8.5 for the years 2050 and 2070 using Support Vector Machine (SVM) algorithm.

Figure 4 .
Figure 4. Global and New Zealand predictions of brown marmorated stink bug habitat suitability under current conditions and emissions trajectories of Representative Concentration Pathways (RCP) 2.6 and RCP 8.5 for the years 2050 and 2070 using Maxent algorithm.

Figure 4 .
Figure 4. Global and New Zealand predictions of brown marmorated stink bug habitat suitability under current conditions and emissions trajectories of Representative Concentration Pathways (RCP) 2.6 and RCP 8.5 for the years 2050 and 2070 using Maxent algorithm.

Figure 5 .
Figure 5. Global and New Zealand predictions of brown marmorated stink bug habitat suitability under current conditions and emissions trajectories of Representative Concentration Pathways (RCP) 2.6 and RCP 8.5 for the years 2050 and 2070 using combined ensemble algorithm (Random Forest, Support Vector Machine and Maxent).

Figure 5 .
Figure 5. Global and New Zealand predictions of brown marmorated stink bug habitat suitability under current conditions and emissions trajectories of Representative Concentration Pathways (RCP) 2.6 and RCP 8.5 for the years 2050 and 2070 using combined ensemble algorithm (Random Forest, Support Vector Machine and Maxent).

Figure 6 .
Figure 6.Centroid shift for Representative Concentration Pathways (RCP) 2.5 and RCP 8.5 from current times to the years 2050 and 2070 using the centroid of reclassified High suitability category areas.

Figure 7 .
Figure 7. Ensembled suitability mapped into the official territories of New Zealand for: current conditions (left); and future climatic conditions (right), for Representative Concentration Pathways (RCP) 8.5 year 2070.

Figure 6 .
Figure 6.Centroid shift for Representative Concentration Pathways (RCP) 2.5 and RCP 8.5 from current times to the years 2050 and 2070 using the centroid of reclassified High suitability category areas.

Figure 6 .
Figure 6.Centroid shift for Representative Concentration Pathways (RCP) 2.5 and RCP 8.5 from current times to the years 2050 and 2070 using the centroid of reclassified High suitability category areas.

Figure 7 .
Figure 7. Ensembled suitability mapped into the official territories of New Zealand for: current conditions (left); and future climatic conditions (right), for Representative Concentration Pathways (RCP) 8.5 year 2070.

Figure 7 .
Figure 7. Ensembled suitability mapped into the official territories of New Zealand for: current conditions (left); and future climatic conditions (right), for Representative Concentration Pathways (RCP) 8.5 year 2070.

Table 1 .
Environmental variables at occurrence points used in the model with final selected variables after checking for cross correlation.

Table 2 .
Evaluation of prediction using Random Forest (RF), Support Vector Machine (SVM) and Maxent models.

Table 2 .
Evaluation of prediction using Random Forest (RF), Support Vector Machine (SVM) and Maxent models.