Biosecurity and post-arrival pathways in New Zealand: relating alien organism detections to tourism indicators

Between-country tourism is established as a facilitator of the spread of invasive alien species; however, little attention has been paid to the question of whether tourism contributes to the arrival and subsequent dispersal of exotic organisms within national borders. To assess the strength of evidence that tourism is a driver for the accidental introducing and dispersal of exotic organisms, we sourced three national databases covering the years 2011 to 2017, namely international and domestic hotel guest nights and national population counts, along with records of exotic organism detections collected by the Ministry for Primary Industries, New Zealand’s government agency that oversees biosecurity. We fitted statistical models to assess the strength of the relationship between monthly exotic organism interception rate, guest nights and population, the latter as a baseline. The analysis showed that levels of incursion detection were significantly related to tourism records reflecting the travel of both international and domestic tourists, even when population was taken into account. There was also a significant positive statistical correlation between the levels of detection of exotic organisms and human population. The core take-home message is that a key indicator of within-country human population movement, namely the number of nights duration spent in specific accommodation, is statistically significantly correlated to the contemporaneous detection of exotic pests. We were unable to distinguish between the effects of international as opposed to domestic tourists. We conclude that this study provides evidence of impact of within-country movement upon the internal spread of exotic species, although important caveats need to be considered.


Introduction
International trade and tourism, while essential to the world's economy, has also been implicated as facilitating the dispersal of exotic species (Early et al. 2016;Seebens et al. 2018). Tourism, in its broadest sense, can provide significant economic gain to a country's GDP, but, if not managed carefully, there are also economic, social, cultural and environmental costs associated with the industry (e.g. Scott et al. 2016;Trivellas et al. 2016;Peeters et al. 2018;Smith et al. 2019). From a biosecurity perspective, the sometimes massive and rapid movement of people associated with international tourism has been implicated in the dispersal of exotic organisms both across and within countries, some of which become invasive (Thuiller et al. 2005;Anderson et al. 2015;Haddaway and Dunn 2015;Early et al. 2016;Hall 2019). Biosecurity failures can have a significant impact on the tourism industry itself, for example, in curtailment of activities once in the country, reducing the value of a country's image to prospective tourists, and a potential reduction in the number of visitors (Blake et al. 2003;Vinson 2013). Exotic species introduced via the tourism pathway can have a direct economic cost, but there are also associated biodiversity losses (Pyšek and Richardson 2010). For countries with a high proportion of endemic biota (such as New Zealand, e.g. McGlone et al. 2001;Lee and Lee 2015), the impacts of exotic species can be significant (Bertram 1999; Barlow et al. 2002;Williams and Timmins 2011).
International tourism has been shown to provide a pathway for the dispersal of many organisms including insects (Russell and Paton 1989;Liebhold et al. 2006), bedbugs (Reinhardt and Siva-Jothy 2007), ticks (Molaei et al. 2019), plant material (Mack and Lonsdale 2001), human diseases (Wilson 1995;Tatem et al. 2006;Khan et al. 2009) and diseased meat products (Pharo 2002). Infested fruit carried by passengers potentially carry unwanted organisms such as fruit fly (Ceratitis and Bactrocera spp.), which could have a significant impact on a country's export fruit industry (SriRamaratnam 1996;Kriticos et al. 2007). Sheridan (1989) found pathogenic fungi on the clothing and baggage of passengers while pockets of clothing have been shown to carry potential risk material including dried and fresh foliage, seeds and feathers (Chirnside et al. 2006). Used tents may carry plant and animal debris, and live insects (Gadgil and Flint 1983). Soiled footwear carried in the luggage of passengers arriving at international airports in New Zealand supported a range of viable bacteria, fungi, seeds and nematodes (McNeill et al. 2011), and included species or strains that were categorized as unwanted organisms under New Zealand's biosecurity regulations. Within a country, contaminated clothing, footwear, camping gear, recreational equipment and vehicles have been shown to provide pathways for dispersing pathogens (Worboys and Gadek 2004;Kidd, et al. 2007), weeds (Whinam et al. 2005;Lloyd et al. 2006;Bouchard et al. 2015) and aquatic organisms (Kilroy and Unwin 2011), into natural and ecologically sensitive environments. Furthermore, the inherent mobility of tourists once within a country's borders (e.g. Forer 2005), also has the potential to facilitate the unintentional transfer of arthropod pests or pathogens from one location to another (Forer and McNeill 2008).
To understand the value of tourism to New Zealand, and thereby associated biosecurity risk, it is worthwhile summarizing some key facts. In the year ended December 2019 there were 3.9 M international visitor arrivals to New Zealand, a 1% increase from the previous year (Stats 2020a). In addition, a further 3.1 M New Zealandresidents returned from overseas holidays (defined as New Zealand residents arriving in New Zealand after an absence of less than 12 months) (Stats 2020b). In the year ended March 2020, total tourism expenditure (both international and domestic) was NZ$42 B, an increase of 2.4% from the previous year, and represented 5.5% of the country's gross domestic product (GDP). For New Zealand, which relies heavily on tourism and primary industries for its economic wealth, biosecurity is strategically important in managing tradeoffs between protecting key economic and environmental assets and encouraging tourism and trade (Parliamentary Commissioner for the Environment 2000).
In this respect, the tension between tourism and biosecurity risk is not unique to New Zealand (e.g. Toral-Granda et al.2017;Melly and Hanrahan 2021), but has been brought into focus with the impact of Covid-19 pandemic and a greater awareness the role tourism plays in dispersing exotic organisms. While the Covid-19 pandemic has had a significant impact on international tourism since the beginning of 2020, the expectation is that in a post-Covid world, there will be a recovery in international travel and renewed growth in global tourism.
In New Zealand, biosecurity monitoring and mitigation of risk at arrival points is a well-established strategy targeting both international and returning New Zealandresident travelers (Jay et al. 2003). But identification and removal of biosecurity risk organisms is not absolute, so tourists (international and returning New Zealanders) may pass through border screening inadvertently carrying undetected risk organisms. Travelers can therefore introduce propagules (sensu Lockwood et al. 2005) in or on their luggage, clothing, and footwear. These propagules (plant pathogens, nematodes, insects, seeds, etc.), can then be deposited at any point along the travel route, depending on activities or events they are undertaking. This could happen when removing boots or jackets from luggage for use while visiting a botanic garden or hiking activity. At this point, a propagule can be deposited at a location where the items were removed from the luggage or along the walking route. The ease for subsequent secondary dispersal may then depend on propagule size and the ability for the exotic organism to be vectored. For example, the introduction of didymo (Didymosphenia geminata, (Lyngbye) M. Schmidt), a freshwater diatom into New Zealand, was strongly linked to anglers arriving from overseas carrying contaminated equipment (Kilroy and Unwin 2011). First found in the southern river systems of the South Island, secondary dispersal was strongly related to human activity, particularly by freshwater anglers (Kilroy and Unwin 2011).
Therefore, understanding the links between international tourist flows once in the country and the potential biosecurity risks that these visitors may present is a new and important area of research. While attempts to visualize tourist movement beyond the port of arrival (either air or sea) within New Zealand, have been made using historical data (e.g. Forer 2005), little is known about subsequent pathways along which international tourists travel. Behavior of tourist flow can differ based on geographical, socioeconomic, demographic, psychographics and behavioral characteristics (Forer 2005;Bigné et al. 2007;Yang et al. 2018). For example, unlike European tourists, Indian and Chinese tourists spend their first few days in Auckland, the main point of arrival into New Zealand before heading elsewhere. In the context of tourism and biosecurity risk, this study sought to address this relationship by using data on (a) New Zealand's monthly hotel guest nights for both international and New Zealand domestic tourists and (b) general population distribution, in relation to biological risk organisms (exotic organisms) detected by New Zealand's biosecurity authority, the Ministry for Primary Industries (MPI). The overall aim was to determine if biosecurity interceptions were best explained by either international or domestic tourist movement within the country, or population density. The broader program would use the results to assist in development of more effective biosecurity risk monitoring and mitigation procedures relating to the different tourist segments. Finally, the information could inform biosecurity authorities on the allocation of resources in relation to other potential pathways (e.g. sea freight).

Materials and methods
We applied a model-comparison approach to assessing the strength of evidence for the competing explanations of the interception patterns. Three data holdings were sourced from the Ministry for Primary Industries (MPI) and Stats NZ Tatauranga Aotearoa (hereafter referred to as Stats NZ). MPI provided the Notification and Investigation Management Application (NIMA) data and Stats NZ, both the monthly hotel domestic and international guest nights data, and annual population data.

Response data (NIMA)
NIMA is the incursion investigation risk identification and reporting framework for notifications to MPI of organisms that may represent a biological risk. The NIMA incursion response data were provided in confidence by MPI and covered the years 2011-2017 (data for earlier years were also provided but not used for the analysis). An incursion is defined by MPI as an exotic organism not previously known to be present in New Zealand, where there is a likelihood that the specimen(s) found is part of a self-sustaining/breeding population. The analysis used the positive records from NIMA as the response variable. A positive also refers to when a risk organism not known to be present in New Zealand is found, but there is no evidence that a selfsustaining / breeding population is present. In this case destroying or treating the risk organism or the risk goods (as the habitat of the organism) removes the threat. The database comprised records of insects, Arachnid spp. (spiders and mites), snails, plants for recording detection areas of exotic organisms b territorial authorities and c region councils from which the annual population datasets were sectioned. New Zealand is divided into 16 regions and 73 territorial authorities. The regions are divided for local government purposes. Territorial authorities are the second tier of local government in New Zealand, below regional councils. Territorial authority districts are not subdivisions of regions, and some of them fall within more than one region. Maps generated using ESRI. ArcGIS Pro. Version 2.7.4. Mar. 6, 2021. https://www.esri.com/en-us/arcgis/products/arcgis-pro/overview  (terrestrial and aquatic), nematodes and microbes (bacteria, fungi and viruses) (all referred hereafter as exotic organisms), their location and date of discovery. Locations were based on the Crosby area codes for recording specimen localities in New Zealand (Crosby et al. 1976). The system comprises 29 geographic areas, with boundaries defined by mountain ranges or rivers, State Highways or straight lines between points (Fig. 1a) (Crosby et al. 1976). A monthly count of NIMA incursion reports is provided in Figure 2, which shows a spike just after the start of each year, corresponding to summer, with smaller winter spikes apparent in some years. NIMA data does not include biological material intercepted at the border such as international airports, seaports or quarantine transitional facilities.

Annual population data
The annual population data were provided by Stats NZ and comprise both city-level and regional annual population data (Fig. 1b, c). In terms of population numbers, the city of Auckland has the highest population, followed by the Canterbury district ( Figure 3).

Tourism data
The tourism data comprised monthly counts of international and domestic visitor nights for accommodation establishments by district for the 2011-2017 period. The accommodation survey collected data on guests (including country of origin) staying in short-term commercial accommodation such as hotels, motels, backpackers, and holiday parks. Domestic data comprises New Zealanders undertaking tourist activities as well as those who may have been away from home for work, family, medical, education and reasons other than simply 'tourism'. Hosted and private accommodation, such as bed and breakfasts and holiday homes, are not included. These include AirBnB, BookaBach, campervans, and friends and family that provide accommodation to both domestic and international guests. While there was no data for 2011 and 2012, this component of accommodation activity was first estimated by Stats NZ in 2013 as 8,4% of the total accommodation industry, rising to 14,5% in 2017 (Grant 2019). Territorial authorities are defined at the meshblock level and represent district and city councils boundaries (Stats 2017). The boundaries of territorial authorities are defined by the 'community of interest', the relevance of the components of the community to each other, and the capacity of the unit to service the community in an efficient manner (Stats 2017).
Auckland was found to dominate domestic occupations, followed by Canterbury (Fig. 4). Domestic tourist nights are much more sharply focused around the New Year than are the international tourist nights, while the pattern for international occupations is much more regular compared to domestic tourist nights. However, international tourist nights show a considerable winter surge in some regions that is not matched by domestic tourists (e.g., Otago Lakes, in particular Queenstown, Figure 4). The total count of international and domestic guest nights are broadly similar across the 20 unique locations (Figure 4).

Managing the spatiotemporal scale
The data represented processes that had been measured at different spatiotemporal levels: the daily (detection), monthly (tourism), and annual (population) levels, and organized variously within the city and district levels. Our goal was to assess whether there is any statistically detectable correlation between the NIMA incursion data and either or both of annual population and monthly tourism data. We chose to construct the model using data corresponding to monthly time-steps, which pick up any seasonal tourism patterns, and at the district or city level.
To complicate matters, the labelling of cities and districts were not consistent within the Stats NZ tourism and population data, respectively. Furthermore, while the district boundaries used in the NIMA incursion data were not the same as applied in the population and tourism databases, there was general alignment with the territorial boundaries used by Stats NZ to segment the latter databases (Fig. 1a, b). We aggregated the datasets to the lowest possible common level of aggregation, leaving 20 distinct locations.
The NIMA data were aggregated to month, and some districts merged to match the population and tourism data. For example, the NIMA data had distinct values for North Canterbury, Mid Canterbury, and South Canterbury, but this level of detail is not supported in the other datasets, so we created a single 'Canterbury' location. The tourism data are reported by month, so no change is needed to the temporal gradient, but as with the other datasets, some merging of district-level data was needed. The population data are annual, so no time changes are needed, and only modest district merges.

Analysis
We applied a forward selection algorithm that starts from a base model and adds (and tests) terms in a curated way. This is because the main alternative, namely backward elimination, involves fitting a complete model and doing so was very time-consuming for these data. The process involved several statistical tests that guided the choices between models. These tests were augmented by other model summary statistics. We compared models using two indices, namely (i) the adjusted R 2 , which can be interpreted as the amount of variation in the response variable that statistically aligns with variation in the predictor variables, adjusted (penalized) to reflect the model size, and (ii) Akaike's Information Criterion (AIC).
The response variable was the number of positive reports each month at a location, which is a non-negative integer. We assumed that the response variable was conditionally Poisson, using a generalized linear modeling approach. We did not consider it safe to assume that the relationship between the candidate predictors and the response variable was a straight line. We fitted a model that allows the relationships to be wiggly, but penalizes the wiggle, so overall it would prefer to be straight, namely an additive model using splines (e.g., Wood 2017). Finally, although the dataset was reasonably large, comprising 1560 monthly observations, it is also highly structured -for example, there are only 20 unique locations (see Figure 1), and the population data are recorded at the year level rather than at the month level. We needed to make the structure of the model match the structure of the data to be confident that the statistical model would discount the data appropriately. We did this by using a mixed-effects model, with year and district random effects. The base model was therefore We applied the following model-fitting approach.
1. We started with a generalized additive mixed-effects model (gamm) that just uses the noted random effects (namely, year and district).
2. We then added annual population as a fixed-effect predictor to account for various levels of otherwise un-measured risk, e.g. sea cargo arrival rates. This term was not formally tested, although its performance will be discussed. This was the base model (above).
3. Next, a penalized smooth function of the sum of domestic and international nights and the difference between domestic and international nights were added, as fixed-effect predictors. Detailed diagnostic checks were made for each model, including spatial and temporal autocorrelation. Checks included: a. Examining a scatterplot of residuals against fitted values to check for obvious lack of fit in the mean or variance model; b. Examining plots of estimated autocorrelation to assess whether the residuals are temporally independent; and c. Adding a smooth surface (a thin-plate spline) in UTM coordinates to see if there is any signal North-South or East-West, which would express in the original models as spatial autocorrelation.
4. This quartet of models (one from step 2 and three from step 3) were then compared, and the comparison interpreted for the statistical information that it provides as to the predictability of positive reports by population and hotel accommodation guest nights.

Results
Model fit statistics are recorded in Table 1. The second row reports the population effect, the third row reports total tourist nights, and the fourth is whether international tourist nights can be distinguished from domestic tourist nights. Both population and total tourist nights are strongly correlated to the number of NIMA incursion reports. Analysis showed that population effect was always monotonically increasing, and either close to linear or linear, and always statistically significant in the model. The total guest nights effect was always linear and increasing and seemed to modify the population effect only very little. Adding international and domestic guest nights improved the Akaike's Information Criterion (AIC). Finally, the difference between international and domestic guest nights effect was flat, suggesting that there is no significant difference in the model based on these data.
The final model of all terms is summarized in Figure 5, which shows that (i) total nights is strongly and linearly related to the natural log of exotic organism reports even when population is considered, and (ii) there is no evidence of any greater risk from international than from domestic nights. A further note on interpretation is that if Table 1. Model fit showing the adjusted R 2 and Akaike's Information Criterion (AIC) in relation to exotic organism interceptions (NIMA reports) in New Zealand. For the AIC values, the lower the number, the better the model fit. The first row reports the base model as defined above; the second is base with (annual) population level added. The third row reports the base model with population and total (monthly) nights of guest nights, and the fourth row includes the previous terms and the difference between international and domestic guest nights. The P-values are generated from the final model in the table (specifically, the full model).  Figure 5. Estimated model effects of the conditional relationship between (i) population and biosecurity incursion reports, (ii) total nights and incursion reports, and (iii) the difference between domestic and international nights and incursion reports. Dashed lines represent approximate 95% confidence limits. the dashed lines intersect (such as in the center and right panel) then the fitted line is straight, otherwise (as in the left panel) it has curvature. The y-axis in each case is labelled with a measure of the magnitude of the curvature needed to capture the relationship. If absent, then the requited curvature is 1, signifying a straight line. Our examination of the model assessment graphics revealed no important caveats (not shown here).

Discussion
The following discussion summarizes the performance of the candidate predictor terms across the set of four nested models that we fitted. There is considerable spatial and temporal variation in the NIMA incursion reports, much of which correlates highly to base human population. The model-fitting exercise shows that there is a clear statistical signal that links reported incursion reports with the hotel guest nights (Table 1). However, it is impossible with the current data to distinguish between the variability that correlates to international as opposed to domestic hotel guest nights. We drew these conclusions using statistical reasoning as follows. Adding the population predictor to the base model greatly enhanced model fit (Table 1). We then added a flexible function of the combined guest nights, that is, the sum of international and domestic guest nights, and found further model improvement, which suggests that variation in guest nights relates to variation in incursion reports that is not otherwise related to population. Finally, we added a flexible function of the subtraction of international from domestic tourist nights. If this term were statistically significant, then this significance would suggest that there is a difference between the effect of domestic as opposed to international tourist nights upon the response variable, and indeed whether the influence of one is greater than the other. No such signal was detected, leading us to conclude that the important apparent relationship is for tourist nights regardless of whether they are international or domestic. The model sketch provided in Figure 5 affirms that both population and total tourist nights are positively correlated to incursion reports, and the model cannot distinguish between the effects of international as opposed to domestic tourist nights. Overall, our analysis showed that tourism, either international or domestic, represents a significant pathway for the introduction and secondary dispersal of biosecurity threats to the extent that this can be established by statistical modeling of an observational study.
We now describe caveats relevant to our interpretation of the model outputs with regards to the underpinning scientific questions. Our goal was to assess the statistical strength of candidate explanatory factors for pest arrival and within-country transport. However, the response variable is the number of exotic organisms detected in the area per month, rather than the number of pests arriving in the area per month. Therefore, we are obliged to assume a tight connection between the arrival of an exotic organism and its detection that amounts to them occurring in the same month. However, this assumption may not always hold; as the research literature shows that a number of historical positives are known to have dispersed undetected, for example emerald ash borer, Agrilus planipennis, (Coleoptera: Buprestidae) in the USA (Siegert et al. 2014) and clover root weevil, Sitona obsoletus Gmelin (Coleoptera: Curculionidae) in New Zealand (Barker et al. 1996). Conversely, early detection of brown marmorated stink bug Halyomorpha halys Stål (Heteroptera: Pentatomidae) in luggage and a hotel room in New Zealand (MPI, unpublished data), exotic fruitfly (Bactrocera and Ceratitis spp.) in surveillance traps (Quilici and Donner 2012), and granulated ambrosia beetle (Xylosandrus crassiusculus, (Motschulsky) (Coleoptera Scolytidae) (Anon 2019), as part of surveillance programs, can improve the probability of determining the pathway for entry as well as improve the probability of eradication if an exotic species were to establish. Positive detection might correlate better with passive surveillance efforts (e.g., citizen science) (Froud et al. 2008;Hester and Cacho 2017) and larger human population centers may have higher probability of detection. We tried to correct for this by including baseline human population in the model, but this assumption could be better tested by a trace-back of each of the reported detections to assess the 'maturity' of the positive at the time of detection. Such an exercise was beyond the remit of our project.
The analysis only considers population count and a measure of within-country tourist activity (monthly number of guest nights). The analysis therefore excludes other potential pathways, including sea freight associated with international trade. The volume of trade imports is generally held to be a more substantial source of biosecurity risk than are international passengers (See Hulme 2009;Sikes et al. 2018). If passenger and freight arrival volumes are correlated, then any potential statistical signal for passengers could be complicated by failing to account for cargo movements. It may be reasonable to believe that arriving freight is correlated with human population concentrations (e.g. Auckland, which is the most populous urban area in New Zealand), however, hitchhiker organisms on freight may have a seasonal pattern that the annual population variable cannot represent. This assumption could be tested by including a candidate predictor that would represent freight activity, for example, monetary value or volumes of imports arriving at both air and sea ports. However, information about subsequent within-country freight movement is not available.
The accommodation survey includes data on short-term commercial accommodation (hotels, motels, backpackers, and holiday parks). Other accommodation types such as 'accommodation-sharing' e.g. AirBnB are not captured, but as noted previously was estimated at c. 8% in 2013 increasing to c. 14% in 2017 of the total accommodation industry (Grant 2019). Therefore, as accommodation-sharing comprised a minor component of the accommodation industry, the biosecurity incursion reports associated with international and domestic visitor nights were highly representative. Within-country tourist accommodation likely comprises a wider variety of activities. For example, some camping grounds and visits to friends and family are not included in the domestic and international guest nights data. Camping trips could make an important difference to the travel statistics, because an unknown proportion will likely involve destinations in more remote or vulnerable areas (see, for example, Runghen at al. 2021). This assumption could be tested by finding further information on campground occupation statistics and recreational vehicle rentals.
The analysis ignores a reasonable supposition that the first few nights for arriving passengers are probably the riskiest from the point of view of the movement of exotic organisms. In the analysis, all nights of accommodation are treated equally. However, the locations of the first few nights for international passengers are likely to be concentrated in areas with high population counts, especially for Auckland, which is New Zealand's main international arrivals airport. On the other hand, analysis of the first seven nights for international passengers shows that they disperse quickly once in the country (Wilson et al. 2018).
The analysis was also unable to discern between New Zealand residents who have arrived from international departure points and New Zealand residents whose travel is purely domestic. However, we consider it reasonable to assume that the influence of returning New Zealanders is relatively negligible in distinguishing between the impact of international and domestic tourism on within-country spread.
The analysis also assumes that the true population data do not change appreciably within the year. Conversely, the other candidate predictors (international and domestic guest nights) both show substantial within-year variation. Therefore, it is possible that the true population data could also change within the year, an assumption that could be assessed if finer-scale data were available.
These results generally support the findings of Edney-Browne et al. (2018), who found that the number of international tourist arrivals to New Zealand, was an important component to explain spatial patterns for establishment of exotic organisms. Conversely, a broader modeling analysis of the major drivers to invasion risk for the "100 among the world's worst invasive alien species", found that socioeconomic variables including human population density, distance to the nearest airport or distance to the nearest seaport, were important contributors to explain the distribution of most taxonomic groups in the list (Bellard et al. 2016).
In conclusion, this analysis using population density and accommodation nights found that the number of reported positive interceptions of exotic organism was significantly positively related to population density and at the same time significantly positively related to total guest nights (combining international and domestic guests). There is no evidence of any difference between international and domestic guests in terms of the relationship with interceptions of exotic organisms. Therefore, we suggest that this study provides conditional evidence that international tourism contributes to the introduction of exotic organism, and within-country movement of both international and domestic tourists aids the secondary dispersal of exotic organisms. While the analyses showed a strong relationship between data for exotic organism interceptions and tourist guest nights, it does not allow us to determine if tourists are also the vector for exotic organisms. However, it may be a reasonable assumption to suggest there is a link which could be investigated. Further research that differentiates the respective role of both tourist segments, and their overall contribution to biosecurity risk in relation to other pathways (e.g. sea freight) for the introduction and dispersal of exotic organisms would also seem warranted. This would contribute to the development of more effective biosecurity risk monitoring and mitigation procedures.
The core take-home message is that anthropogenic movements associated with tourism correlate with detection of exotic organisms in New Zealand. The results also reinforce the need for biosecurity authorities to continue to allocate resources to managing the tourism pathway.