Spatial Distribution of the Mexican Daisy, Erigeron karvinskianus , in New Zealand under Climate Change

: The invasive species Erigeron karvinskianus or Mexican daisy is considered a signiﬁcant weed that impacts native forest restoration efforts in New Zealand. Mapping the potential distribution of this species under current and future predicted climatic conditions provides managers with relevant information for developing appropriate management strategies. Using occurrences available from global and local databases, spatial distribution characteristics were analyzed using geostatistical tools in ArcMap to characterize current distribution. Species distribution modeling (SDM) using Maxent was conducted to determine the potential spatial distribution of E. karvinskianus worldwide and in New Zealand with projections into future climate conditions. Potential habitat suitability under future climatic conditions were simulated using greenhouse gas emission trajectories under the Representative Concentration Pathway (RCP) models RCP2.6, RCP4.5, RCP6.0 and RCP8.5 for years 2050 and 2070. Occurrence data were processed to minimize redundancy and spatial autocorrelation; non-correlated environmental variables were determined to minimize bias and ensure robust models. Kernel density, hotspot and cluster analysis of outliers show that populated areas of Auckland, Wellington and Christchurch have signiﬁcantly greater concentrations of E. karvinskianus . Species distribution modeling results ﬁnd an increase in the expansion of range with higher RCP values, and plots of centroids show a southward movement of predicted range for the species.


Introduction
Mexican Daisy (Erigeron karvinskianus DC.) is reported as an invasive plant species throughout warm-temperate, subtropical and tropical parts of the world [1], including India [2], Africa [3], Australia [4], Japan [5], northwest Himalayas [6] and Europe [7,8]. Classified as an invasive weed, it is an identified problem in Hawaii [9,10], La Réunion [11] and New Zealand [12,13]. The species was first recorded in New Zealand in 1940 [12,13] and is now found to be established and thriving in North Island, South Island and Stewart Island [13,14], with numerous patches found in the relatively populated Auckland Region [12,15].
The species is known to be associated with anthropogenic disturbance, particularly in habitats such as rock walls, archaeological sites, roadside banks, and wasteland [1,4,8,13,[16][17][18] but is also increasingly invading native ecosystems including native and replanted forests, wetlands and riparian zones [11,16,19]. It can form dense mats that can rapidly stifle young plants [1,9], with serious impacts on projects involving revegetation or replanting of areas. This in turn may require significant resources for clearing, preparation and maintenance to ensure the growth and survival of replanted species. The plant's ability to spread at a rapid rate when introduced further increases its importance in

Materials and Methods
Occurrence data were sourced from the Global Biodiversity Information Facility (GBIF; http: //data.gbif.org). As a database with data sourced from a wide variety of museums, projects, networks and herbaria worldwide with varying collection quality, there are recognized challenges in its use for SDM [40]. However, the ability to access and download a set of global georeferenced species occurrences makes it a valuable resource, particularly for species with no other alternative source of information. The downloaded distribution data were checked for consistency, and redundant, unnecessary and non-georeferenced data were removed. The processed Microsoft Excel data were imported into ArcMap, and the resulting distribution map was used for checking obvious errors such as redundancies and data at improbable locations, such as bodies of water or extreme latitudes.
Using the Kernel Density tool of ArcMap, a density map was created to show areas where there are greater concentrations of reported occurrences. This tool calculates magnitude of occurrences per unit area and fits a suitably smooth surface over the extent of input data points using the distances between each location. The raster output surface aids the visualization of the distribution by highlighting areas that have much greater or much lower concentration of occurrences. A planar method for generating the surface was the option used in the calculations [41,42].
The Getis-Ord Gi* [43,44] statistic implemented in the Optimized Hotspot tool of ArcMap determined characteristics of occurrences in terms of statistically significant intensities of occurrence points. The Hotspot tool aggregates occurrences into uniformly-sized cells and calculates the z-statistic with an accompanying p-value that indicates whether the cells are hotspots, coldspots or have no significance. Positive z-statistics identify cells that have high numbers of occurrences, and designates them as hot spots when the p-value is significant [45]. At the opposite scale, cells with negative z-statistics represent low numbers of occurrences and are designated as cold spots if the calculated p-value is significant [45]. To characterize the neighborhood and identify areas with characteristics of statistically significant outliers in the occurrence cells, the Optimized Cluster and Outlier tool of ArcMap that implements Anselin's Local Moran's I [46] was used. This tool produces cells that aggregate occurrences and use counts per cell as the variable to determine a z-score to depict outlier significance. High or low values of z-scores result in cells depicted as High-High or Low-Low cells in the map. Other results include the Low-High and High-Low cells showing the nature of the surrounding cells in relation to each other. Cells not determined by the tool to have statistically significant outliers are characterized as not significant [47].
Species distribution modeling requires a set of occurrence data and environmental variables as input. The tool Maxent (v3.4.1), implementing the maximum entropy machine learning approach, was used. Maxent is a popular tool finding widespread usage because of its performance compared to other approaches [48,49]. Maxent is based on maximum entropy that calculates inferences or predictions even under conditions where the information is not complete [50][51][52]. In SDM applications, the maximum entropy approach estimates a target distribution by calculating the probability distribution that is nearest to a uniform one under constraints determined from available data [50]. Significant applications of Maxent for modeling plant invasive species include modeling of invasive trees worldwide [53], invasive weeds in Australia [30], the Chinese fan palm [39], and the tallow tree in the Himalayas [54].
Environmental data consisted of rasters available from the WorldClim (v1.4) database that consist of 19 Bioclim environmental variables (http://worldclim.org) [55] commonly used for SDM in a wide variety of applications [56][57][58]. Current Bioclim rasters (averaged between 1960 and 1990) and available future climatic data were downloaded for 2050 and 2070 at 30-second resolution (http: //www.worldclim.org/CMIP5v1). These future climatic data describe different levels of greenhouse gas trajectories as defined in the 5th IPCC report. The downloaded raster files represent Representative Concentration Pathways (RCP) with values of 2.6, 4.5, 6.0, and 8.5 (higher numbers represent greater greenhouse gas emissions, with RCP2.6 resulting in global warming by less than 2 degrees and RCP8.5 likely to lead to 4 degrees warming by 2100) for the years 2050 and 2070 [59]. The downloaded occurrence data were rarefied using the SDMToolbox v2.2 [60] in ArcMap with multiple presence points within specified range distances reduced to a single representative occurrence in order to minimize spatial autocorrelation which can affect model robustness [44,61,62]. Five sets of occurrence data were generated consisting of rarefied sets at distances of 1, 5, 10, 20, 50, and 100 km between samples. Each rarefied set was modelled in preliminary Maxent runs to determine which data set to use based on its AUC [63] performance for both worldwide and New Zealand occurrence data sets. The higher AUC score was used to select the rarefied data set for the global model, whereas for the New Zealand model, the minimal difference in AUC between the 1 km and 5 km resolutions resulted in the decision to use the 5 km data with a smaller number of samples to further minimize spatial autocorrelation [39]. The sample sizes for both models are also well above the minimum of 30 (Table 1) recommended for a range of SDM algorithms [64]. Maxent outputs a suitability map by modeling the distribution of a specific species with a selected number of environmental variables. The worldwide model was used to determine the non-correlated variables to be used for the New Zealand model. The New Zealand current model in turn was projected into future environmental conditions using RCP bioclimatic rasters. The Minimum Training Presence (MTP) threshold rule was applied in order to produce a binary raster representing the presence/absence or range of a species [65]. The MTP as a threshold in Maxent was used [66][67][68][69] primarily because of its representation of the simple interpretation of an ecological relationship between areas that are at least suitable to locations where presence is recorded. Results of the preliminary Maxent run identified the 100 k data set as most suitable for the global model and the 5 k data set for the New Zealand model (Table 1). Bioclim rasters were further processed as it was expected that many of the climatic variables were highly correlated [60], making it difficult to determine the effects of individual variables on the Maxent model. These consisted of 19 parameters in Bioclim clipped to the New Zealand area and rasters representing elevation (https://data.linz.govt.nz/layer/51768-nz-8m-digital-elevation-model-2012/) and land cover (https://earthexplorer.usgs.gov/) [70,71], using the Biosphere 2 classification of 11 landcover types (1 Broadleaf Evergreen Trees; 2 Broadleaf Deciduous Trees; 3 Broadleaf and Needleleaf Trees; 4 Needleleaf Evergreen Trees; 5 Needleleaf Deciduous Trees; 6 Short Vegetation/Grassland; 7 Shrubs with Bare Soil; 8 Dwarf Trees and Shrubs; 9 Agriculture/Grassland; 10 Water, Wetlands; 11 Ice/Snow) and human population (https://koordinates.com/from/datafinder.stats.govt.nz/layer/ 8437/). Landcover, elevation and human population variables were included in the model together with bioclimatic factors to determine their effects on the Maxent jackknife outputs. Elevation was selected primarily because of observed range shifts towards higher elevation in climate changes studies [72] and to show results of variable effects on the models in the jackknife of the Maxent model while considering that New Zealand's mountain areas contain significant habitats that may be impacted. The values of the environmental variables at each occurrence point of the selected rarefied data set were determined using the ArcMap tool Extract Multi-values from rasters ( Table 2). We used SDMtoolboxv2.2 in ArcMap to check for cross-correlation of all the environmental variables at all occurrence points and determine which environmental variables to use based on the correlation coefficient [60]. Variables not correlated at coefficient less than or equal to 0.7 were included and further checked to include only those with a variance inflation factor (VIF) of less than 10 [73]. This is consistent with results that show that using simpler models with a small number of variables improves the results of projecting species distribution models [58,[74][75][76]. Cross correlation of the Bioclim data resulted in the 4 layers plus landcover, population and elevation selected for use in both the global and New Zealand models ( Table 2). The global and New Zealand species distribution models were run with the 100 k and 5 k rarefied data sets, respectively. Environmental variables consisted of non-correlated current Bioclim environmental layers BIO1 (Annual Mean Temperature), BIO2 Mean Diurnal Range (Mean of monthly (max temperature-min temp)), BIO12 (Annual Precipitation), and BIO15 (Precipitation Seasonality (Coefficient of Variation), elevation, land cover and human population. The New Zealand model was projected into future climatic conditions using the downloaded RCP trajectories for the years 2050 and 2070.
Using raster calculations, subtracting the binary thresholded future rasters from current rasters resulted in four categories: Expansion in range (absence in current, presence in future), no occupancy Climate 2019, 7, 24 6 of 20 or (absence in both current and future), occupancy or continuing (presence in current and future) and contraction in range (presence in current, absence in future). These resulting range maps provide valuable insight into areas of the country where change in range occurs. The maps also show effects of different RCP values on the range expansion or contraction of the species. The direction of presence area centroid movement was determined in ArcMap, again using the tool Distribution Changes Between Binary SDMs in the SDM toolbox. Centroid movements are plotted as vectors and represent the magnitude and direction of the movement of predicted distribution ranges. Each centroid shift vector was produced by pairwise subtraction between binary current and future rasters and also between the 2050 and 2070 binary rasters. The speed of the movement of the centroid is then calculated from the length of the centroid shift and the time difference between each pair of presence/absence rasters.

Results
Global spatial distribution characterization showed most of the occurrences were found in the southern region of Mexico, a large area of Europe, and the east coasts of Australia and New Zealand ( Figure 1). The Optimized Hotspot Tool in ArcMap depicted the native range of the plant in Mexico and New Zealand as significant global hotspots. The rest of the world did not exhibit any significant hotspots in the generated cells analyzed. In terms of the cluster and outlier characteristics of the global distribution, the native range of the plant in Mexico showed several significant High-High cluster cells as well as Low-High outlier characterized cells. New Zealand shows several Low-High outlier cells, whereas a few High-Low outlier cells are found scattered over the world ( Figure 2). For the New Zealand distribution, the downloaded occurrence data plus the survey in the Auckland Region show that the Mexican daisy is currently prevalent in regions surrounding Auckland, Wellington and Christchurch and scattered throughout the length of North Island, and South Island, particularly at the northern end and in the Christchurch region ( Figure 2). This is confirmed by kernel density modeling that shows Auckland, Wellington, the top of South Island and Christchurch to have higher densities. Areas at the eastern tip of North Island and the southwestern area of South Island do not have reported occurrences and are excluded in the resulting kernel density raster ( Figure 3). Hotspot analysis using the Getis-Ord Gi* statistic show that the center of Auckland has statistically significant hotspots. Cluster and outlier analysis find Auckland with several significant Low-High outliers. Several spots of High-Low outliers are found throughout North Island and only one High-Low outlier cell in South Island was found. (Figure 3)  Limiting the occurrence data set to the Auckland Region alone, there is a high density of occurrences seen particularly around the western and coastal areas. When the Optimized Hotspot For the New Zealand distribution, the downloaded occurrence data plus the survey in the Auckland Region show that the Mexican daisy is currently prevalent in regions surrounding Auckland, Wellington and Christchurch and scattered throughout the length of North Island, and South Island, particularly at the northern end and in the Christchurch region ( Figure 2). This is confirmed by kernel density modeling that shows Auckland, Wellington, the top of South Island and Christchurch to have higher densities. Areas at the eastern tip of North Island and the southwestern area of South Island do not have reported occurrences and are excluded in the resulting kernel density raster ( Figure 3). Hotspot analysis using the Getis-Ord Gi* statistic show that the center of Auckland has statistically Climate 2019, 7, 24 9 of 20 significant hotspots. Cluster and outlier analysis find Auckland with several significant Low-High outliers. Several spots of High-Low outliers are found throughout North Island and only one High-Low outlier cell in South Island was found ( Figure 3).
Climate 2018, 6, x FOR PEER REVIEW 9 of 20 tool was run, no significant hotspots were found in the distribution. Cluster and Outlier analysis, however, finds High-Low clusters in the northern and western areas of the city (Figure 3). Species distribution modeling results from global Maxent modeling with the selected data set show the native range in Mexico and most reported invasive range of the species to be the most suitable under current conditions. Further, a significant area of New Zealand was found to be highly suitable. Coastal areas of Europe, the southern coast of Australia, the Japanese archipelago, central China, the eastern United States and southeastern Central America are also highly suitable. The presence/absence map that resulted with a threshold set at MTP level shows that a significant area of the world is predicted to be favorable for its presence or can be considered as the range of the species under current environmental conditions ( Figure 4). Limiting the occurrence data set to the Auckland Region alone, there is a high density of occurrences seen particularly around the western and coastal areas. When the Optimized Hotspot tool was run, no significant hotspots were found in the distribution. Cluster and Outlier analysis, however, finds High-Low clusters in the northern and western areas of the city (Figure 3).
Species distribution modeling results from global Maxent modeling with the selected data set show the native range in Mexico and most reported invasive range of the species to be the most suitable under current conditions. Further, a significant area of New Zealand was found to be highly suitable. Coastal areas of Europe, the southern coast of Australia, the Japanese archipelago, central China, the eastern United States and southeastern Central America are also highly suitable. The presence/absence map that resulted with a threshold set at MTP level shows that a significant area of the world is predicted to be favorable for its presence or can be considered as the range of the species under current environmental conditions (Figure 4). Climate 2018, 6, x FOR PEER REVIEW 10 of 20 For the New Zealand Maxent model, jackknife analysis of variable contributions shows that human population has the greatest percentage contribution followed by land cover and BIO1 (average annual temperature). However, when permutation importance is considered, BIO1 is more important, followed by land cover and elevation. Response curves provide an idea of the variable value range or category that results in the greatest model response. The response curve for BIO1 (annual average temperature range) shows that the model has a higher response at the higher values of annual average temperature. Middle values of BIO15 (Precipitation seasonality) and lower values of BIO2 (Mean diurnal temperature range) also have more influence in the model. Lower elevations have greater effect compared to higher elevations, whereas BIO12 (average annual rainfall) has the greatest effect at lower averages, particularly when it is the only corresponding variable. For land cover, categories 1 (broadleaf evergreen forest) and 7 (shrubs with bare soil) consistently show greater effect in both response graph types (Table 2).  For the New Zealand Maxent model, jackknife analysis of variable contributions shows that human population has the greatest percentage contribution followed by land cover and BIO1 (average annual temperature). However, when permutation importance is considered, BIO1 is more important, followed by land cover and elevation. Response curves provide an idea of the variable value range or category that results in the greatest model response. The response curve for BIO1 (annual average temperature range) shows that the model has a higher response at the higher values of annual average temperature. Middle values of BIO15 (Precipitation seasonality) and lower values of BIO2 (Mean diurnal temperature range) also have more influence in the model. Lower elevations have greater effect compared to higher elevations, whereas BIO12 (average annual rainfall) has the greatest effect at lower averages, particularly when it is the only corresponding variable. For land cover, categories 1 (broadleaf evergreen forest) and 7 (shrubs with bare soil) consistently show greater effect in both response graph types (Table 3). Using the 5km rarefied data for New Zealand species distribution modeling reflects the distribution pattern (hotspots and cluster/outlier analysis) that resulted from the geostatistical tools in ArcMap. Higher suitability for the Auckland, Wellington and Christchurch regions is evident in the maps. When projected into future trajectories, similar distributions are found with some changes in the distribution ( Figure 5). Using the 5km rarefied data for New Zealand species distribution modeling reflects the distribution pattern (hotspots and cluster/outlier analysis) that resulted from the geostatistical tools in ArcMap. Higher suitability for the Auckland, Wellington and Christchurch regions is evident in the maps. When projected into future trajectories, similar distributions are found with some changes in the distribution ( Figure 5). Using the 5km rarefied data for New Zealand species distribution modeling reflects the distribution pattern (hotspots and cluster/outlier analysis) that resulted from the geostatistical tools in ArcMap. Higher suitability for the Auckland, Wellington and Christchurch regions is evident in the maps. When projected into future trajectories, similar distributions are found with some changes in the distribution ( Figure 5). Using the 5km rarefied data for New Zealand species distribution modeling reflects the distribution pattern (hotspots and cluster/outlier analysis) that resulted from the geostatistical tools in ArcMap. Higher suitability for the Auckland, Wellington and Christchurch regions is evident in the maps. When projected into future trajectories, similar distributions are found with some changes in the distribution ( Figure 5). BIO1 17 . Using the 5km rarefied data for New Zealand species distribution modeling reflects the distribution pattern (hotspots and cluster/outlier analysis) that resulted from the geostatistical tools in ArcMap. Higher suitability for the Auckland, Wellington and Christchurch regions is evident in the maps. When projected into future trajectories, similar distributions are found with some changes in the distribution ( Figure 5). Using the 5km rarefied data for New Zealand species distribution modeling reflects the distribution pattern (hotspots and cluster/outlier analysis) that resulted from the geostatistical tools in ArcMap. Higher suitability for the Auckland, Wellington and Christchurch regions is evident in the maps. When projected into future trajectories, similar distributions are found with some changes in the distribution ( Figure 5). Using the 5km rarefied data for New Zealand species distribution modeling reflects the distribution pattern (hotspots and cluster/outlier analysis) that resulted from the geostatistical tools in ArcMap. Higher suitability for the Auckland, Wellington and Christchurch regions is evident in the maps. When projected into future trajectories, similar distributions are found with some changes in the distribution ( Figure 5). Using the 5km rarefied data for New Zealand species distribution modeling reflects the distribution pattern (hotspots and cluster/outlier analysis) that resulted from the geostatistical tools in ArcMap. Higher suitability for the Auckland, Wellington and Christchurch regions is evident in the maps. When projected into future trajectories, similar distributions are found with some changes in the distribution ( Figure 5). Using the 5km rarefied data for New Zealand species distribution modeling reflects the distribution pattern (hotspots and cluster/outlier analysis) that resulted from the geostatistical tools in ArcMap. Higher suitability for the Auckland, Wellington and Christchurch regions is evident in the maps. When projected into future trajectories, similar distributions are found with some changes in the distribution ( Figure 5). Using the 5km rarefied data for New Zealand species distribution modeling reflects the distribution pattern (hotspots and cluster/outlier analysis) that resulted from the geostatistical tools in ArcMap. Higher suitability for the Auckland, Wellington and Christchurch regions is evident in the maps. When projected into future trajectories, similar distributions are found with some changes in the distribution ( Figure 5).  Using the 5km rarefied data for New Zealand species distribution modeling reflects the distribution pattern (hotspots and cluster/outlier analysis) that resulted from the geostatistical tools in ArcMap. Higher suitability for the Auckland, Wellington and Christchurch regions is evident in the maps. When projected into future trajectories, similar distributions are found with some changes in the distribution ( Figure 5). Using the 5km rarefied data for New Zealand species distribution modeling reflects the distribution pattern (hotspots and cluster/outlier analysis) that resulted from the geostatistical tools in ArcMap. Higher suitability for the Auckland, Wellington and Christchurch regions is evident in the maps. When projected into future trajectories, similar distributions are found with some changes in the distribution ( Figure 5). Using the 5km rarefied data for New Zealand species distribution modeling reflects the distribution pattern (hotspots and cluster/outlier analysis) that resulted from the geostatistical tools in ArcMap. Higher suitability for the Auckland, Wellington and Christchurch regions is evident in the maps. When projected into future trajectories, similar distributions are found with some changes in the distribution ( Figure 5). Using the 5km rarefied data for New Zealand species distribution modeling reflects the distribution pattern (hotspots and cluster/outlier analysis) that resulted from the geostatistical tools in ArcMap. Higher suitability for the Auckland, Wellington and Christchurch regions is evident in the maps. When projected into future trajectories, similar distributions are found with some changes in the distribution ( Figure 5).
Using the 5 km rarefied data for New Zealand species distribution modeling reflects the distribution pattern (hotspots and cluster/outlier analysis) that resulted from the geostatistical tools in ArcMap. Higher suitability for the Auckland, Wellington and Christchurch regions is evident in the maps. When projected into future trajectories, similar distributions are found with some changes in the distribution ( Figure 5). Thresholded rasters show small differences between the current model and the future distribution ( Figure 6). A noticeable change in the predicted presence areas is found in the higher emission RCP8.5 for 2070 where large areas of North Island show a greater absence prediction compared to the rest of the thresholded maps. The RCP4.5 and RCP6.0 maps show greater predicted absence areas compared to the current habitat suitability map.  Thresholded rasters show small differences between the current model and the future distribution ( Figure 6). A noticeable change in the predicted presence areas is found in the higher emission RCP8.5 for 2070 where large areas of North Island show a greater absence prediction compared to the rest of the thresholded maps. The RCP4.5 and RCP6.0 maps show greater predicted absence areas compared to the current habitat suitability map. Thresholded rasters show small differences between the current model and the future distribution ( Figure 6). A noticeable change in the predicted presence areas is found in the higher emission RCP8.5 for 2070 where large areas of North Island show a greater absence prediction compared to the rest of the thresholded maps. The RCP4.5 and RCP6.0 maps show greater predicted absence areas compared to the current habitat suitability map.  When the change in range characteristics of the thresholded presence/absence maps is calculated between current conditions and the RCP rasters for 2070, the contraction of the areas of species range is less than the area expansion, except for RCP8.5, where the opposite is true. The increase in the expansion range is shown to be directly proportional with increasing RCP values, with RCP8. 5 showing an expansion in range almost double that of RCP2.6 (22,749 km 2 vs. 12,089 km 2 ). On the other hand, the areas that contract in RCP8.5 are much greater than the other scenarios combined. There is also a decrease in areas that show no change or presence in both rasters with the higher RCP values (Figure 7 and Table S1).
When the change in range characteristics of the thresholded presence/absence maps is calculated between current conditions and the RCP rasters for 2070, the contraction of the areas of species range is less than the area expansion, except for RCP8.5, where the opposite is true. The increase in the expansion range is shown to be directly proportional with increasing RCP values, with RCP8. 5 showing an expansion in range almost double that of RCP2.6 (22,749km 2 vs. 12,089km 2 ). On the other hand, the areas that contract in RCP8.5 are much greater than the other scenarios combined. There is also a decrease in areas that show no change or presence in both rasters with the higher RCP values (Figure 7 and Table S1). The range maps resulting from subtracting the thresholded future raster from the current raster show that North Island becomes less favorable for the Mexican daisy by the year 2070 for all RCP trajectories. An obvious range expansion in South Island is evident with increasing RCP values [77]. When centroid shifts were determined between the current and future as well as from 2050 to 2070 thresholded rasters, all shifts showed a movement and direction to the southwest, mainly along the length of New Zealand's longitudinal axis (Figure 8). The magnitude of the centroid shift is also greater with higher RCP values, with RCP8.5 showing the highest magnitude. Centroid shifts for RCP4.5 is greater than RCP2.6, with the year-current-to-2070 shift greater than the shift from current to 2050. (Figure 8). When the range shift was calculated using the movement of the centroid, the values in were 0.57km/yr (RCP2.6), 0.85km/yr (RCP4.5), 1.06km/yr (RCP6.0) and 1.41km/yr (RCP8.5). The range of this movement is greater than the post-glacial migration rate of 140 plants estimated at values from 0.01 to 0.44 km/yr, with a mean of 0.12 km/yr [78]. The range maps resulting from subtracting the thresholded future raster from the current raster show that North Island becomes less favorable for the Mexican daisy by the year 2070 for all RCP trajectories. An obvious range expansion in South Island is evident with increasing RCP values [77]. When centroid shifts were determined between the current and future as well as from 2050 to 2070 thresholded rasters, all shifts showed a movement and direction to the southwest, mainly along the length of New Zealand's longitudinal axis (Figure 8). The magnitude of the centroid shift is also greater with higher RCP values, with RCP8.5 showing the highest magnitude. Centroid shifts for RCP4.5 is greater than RCP2.6, with the year-current-to-2070 shift greater than the shift from current to 2050.

Discussion
This work produced maps depicting spatial distribution and potential suitability of the Mexican daisy at global and New Zealand scales. The global maps show the widespread distribution of the Mexican daisy at its native and invaded areas. Regions with similar temperatures to its native range and coastal areas consistently show higher levels of suitability for the species. At a global scale, New Zealand shows very high levels of suitability and predicted presence/absence for the species. The species appears to be most prevalent in coastal regions; in addition to most of New Zealand it can also be seen along the southern coast of Australia, South America, Japan, and some areas in Europe.
There are a few outliers where the species has been noted growing further inland, such as Africa and China. Outside of its native range, Optimized Hot Spot Analysis as well as CLOA show that New Zealand is the only country in the world with significant hotspots and significant outliers, a result that may be due to the number of occurrences reported in GBIF and that could indicate an inherent spatial bias due to uneven reporting [79]. This is also consistent with the country's top 2 ranking in the number of invasive species recorded [80].
Recognized qualifications relevant in the use of SDM [81,82 ] related to this effort include the nature of occurrence data sourced from GBIF, use of a single performance metric in model evaluation and model parameters that may need iterative values to produce better results. In the absence of any other source of georeferenced location data, the rarefication and data processing conducted may also be considered as a form of subsampling that preserves the extent and contribution towards better predictive model [79] and removes some of the effects of sampling bias [83] because the source data is sufficiently large [81]. Further work to improve this aspect of modeling is certainly warranted. Using other evaluation metrics such as True Skills Statistics (TSS) and kappa [39,62] should increase confidence in the selection of models. For the model parameters specifically on threshold choice, related studies testing different threshold values derived from several model parameters such as sensitivity and specificity show more superior threshold options than the one used in this study [84,85]. Other threshold options, such as the maximum sum of sensitivity and specificity, that were shown to perform better [86,87] could therefore be tested and compared for use in future iterations of the species distribution modeling for the Mexican daisy.

Discussion
This work produced maps depicting spatial distribution and potential suitability of the Mexican daisy at global and New Zealand scales. The global maps show the widespread distribution of the Mexican daisy at its native and invaded areas. Regions with similar temperatures to its native range and coastal areas consistently show higher levels of suitability for the species. At a global scale, New Zealand shows very high levels of suitability and predicted presence/absence for the species. The species appears to be most prevalent in coastal regions; in addition to most of New Zealand it can also be seen along the southern coast of Australia, South America, Japan, and some areas in Europe.
There are a few outliers where the species has been noted growing further inland, such as Africa and China. Outside of its native range, Optimized Hot Spot Analysis as well as CLOA show that New Zealand is the only country in the world with significant hotspots and significant outliers, a result that may be due to the number of occurrences reported in GBIF and that could indicate an inherent spatial bias due to uneven reporting [79]. This is also consistent with the country's top 2 ranking in the number of invasive species recorded [80].
Recognized qualifications relevant in the use of SDM [81,82] related to this effort include the nature of occurrence data sourced from GBIF, use of a single performance metric in model evaluation and model parameters that may need iterative values to produce better results. In the absence of any other source of georeferenced location data, the rarefication and data processing conducted may also be considered as a form of subsampling that preserves the extent and contribution towards better predictive model [79] and removes some of the effects of sampling bias [83] because the source data is sufficiently large [81]. Further work to improve this aspect of modeling is certainly warranted. Using other evaluation metrics such as True Skills Statistics (TSS) and kappa [39,62] should increase confidence in the selection of models. For the model parameters specifically on threshold choice, related studies testing different threshold values derived from several model parameters such as sensitivity and specificity show more superior threshold options than the one used in this study [84,85]. Other threshold options, such as the maximum sum of sensitivity and specificity, that were shown to perform better [86,87] could therefore be tested and compared for use in future iterations of the species distribution modeling for the Mexican daisy.
Looking at the relationship between the distribution of the species at both global and New Zealand scales can assist in managing the spread of the species within New Zealand. Currently the Mexican daisy is growing in many major regions across the globe and New Zealand, most notably around the major urban centers Auckland, Wellington, Christchurch and coastal regions of South Island. This is consistent with the characteristic of the Mexican daisy in forest areas to be associated with the size and proximity of settlements [18]. The modeling results provide information for conservation and management planning for areas of considerable importance most likely to be impacted by the spread and establishment of the species. The information inherent in the maps showing the geographical distribution and results of the geostatistical processing in ArcMap as well as SDM provides graphical knowledge contributing to the effectiveness of related control measures by providing a guide to areas for focusing resources or prioritizing mitigation or eradication measures that are required [88]. Regions with the higher current density show changes in suitability under future climatic conditions. Further exploration of other environmental and even socio-economic characteristics in areas such as Auckland, Wellington and South Island will refine the model and contribute to an increased understanding of the species' ability to establish, grow and spread successfully in these regions. The spread into southern regions in New Zealand is of concern especially in novel areas where native plants are at risk from the competition of the species invading its ecological niche and other associated impacts on the conservation of biodiversity for New Zealand in general. The decrease in habitat suitability overall is consistent with modeling work on invasive species in Australia that compares current and future conditions of hundreds of invasive plants [89].
There is evidence that the species will increase its range in South Island, as well as regions in the south of North Island. This is consistent with the report in the Himalayas that include the Mexican daisy among 11 invasive species modeled similarly and depicting that such will spread to higher elevation and latitudes with global warming [39]. Species distribution modeling predicts that the Mexican Daisy will spread to other areas where it may not currently be a concern under conditions of increasing greenhouse gas emissions.
This information will allow managers of those areas becoming more suitable, particularly in the southern areas of the country, to prioritize or increase the risk rating of the Mexican daisy for their respective regions. Conversely, in areas where the model shows a contraction of range, less focus on the species would reduce its risk rating, providing information for allocating resources to other more risky plants or organisms [90].

Conclusions
The results of this work provide an overview of the potential impacts and spread of the Mexican daisy in New Zealand. The combination of geostatistical processing and SDM provided several predictive maps that can be used to assesses measures for controlling the species. The resulting maps provided clear indications where the species will most likely spread over the next several years, primarily because of changing and warming climate regime as represented in the environmental variables used. This provides a basis and starting point from which to formulate a management plan and also point out for the need for additional research work, including enhancing the data needed for more robust and relevant modeling.