A distribution database of terrestrial isopods was developed in Japan that includes distribution data of 17,412 specimens based on 159 publications, a database of specimens deposited in Japanese museums (S-Net; National Museum of Nature and Science 2014) and personal collections of the first author. The distribution data of the seven species are available in Suppl. material 1 and all the distribution data is visualised by the web system developed for this study (Suppl. material 2). In addition, occurrence data of the seven species were obtained for their native region (Europe) and another invaded region (North America) from the Global Biodiversity Information Facility (GBIF; GBIF.org 2017) because sufficient numbers of occurrence data were reported from both regions. The area of Europe is defined as N35°–N70° and W11°–E23° and the North American region was defined as N10°–N55° and W50°–W135°. The original regions of the seven species examined were unclear, but A. nasatum may have originally been distributed in Italy, southern France and Spain and A. vulgare and P. pruinosus may be native to the Mediterranean region (Cochard et al. 2010). For the present study, the authors tentatively considered the European continent as the original distribution area of these taxa, although some species were likely foreign in some parts of Europe (Cochard et al. 2010). In addition, P. scaber was described as a subspecies based on specimens collected from Hokkaido, northern Japan (Verhoeff 1928), but this taxonomic classification was dubious. For example, genetic characters of two specimens from Japan were almost identical to those of specimens from Europe and the United States (Karasawa et al. 2017); this species has been treated as exotic in Japan (Nunomura 2007). Duplicate locations were excluded for further analyses.

Bioclimatic variables (19) were used from the WorldClim version 2.0 dataset at five arc-minute resolution for these analyses (Fick and Hijmans 2017). These variables included temperature and precipitation data, which have been recommended for and widely applied to species distribution modelling (e.g. Broennimann et al. 2007, Kumar et al. 2015, Zhu et al. 2017). Highly collinear variables (Pearson correlation coefficient, |r| ≥ 0.75) were removed and only one variable from a set of highly correlated variables was included in each model (Suppl. materials 3, 4).

The environmental niche models (ENMs) were constructed using Maxent version 3.4.1 (Phillips et al. 2017), which employs the principle of maximum entropy to discriminate occurrence records of species to environmental variables to estimate a species’ potential geographic distribution (Phillips et al. 2006). The Maxent programme has the following advantages compared with other methods: the technique can be performed using presence data only (Phillips et al. 2006) and is relatively robust with small sample sizes (Hernandez et al. 2006; Wisz et al. 2008). Two models were constructed: 1) the regional model developed based on occurrence locations in the invasive region (Japan) and 2) the global model developed based on all occurrence locations in Japan, Europe and North America. The occurrence data may hold sampling bias because they were not collected randomly and were obtained from several sources. To enhance the performance for constructing ENMs, background data were created with the same underlying bias for the occurrence data (target group background; Philips et al. 2009). A bias surface was generated using a kernel density estimation using the sm package (Bowman and Azzalini 2014, the script was according to Fitzpatrick et al. 2013) in R version 3.4.3 (R Development Core Team 2017). Then, 10,000 background points were randomly selected as weighed by the sampling bias surface. Linear and quadratic features were used as functional forms to describe relationship probability of occurrence and environmental data, to prevent over-fitting of the training data (Syfert et al. 2013, 2014). Model accuracy was evaluated based on the area under the curve (AUC). AUC values indicated whether model predictions were better than random; values < 0.5 denote worse than random performance; 0.5–0.7 indicate poor performance; 0.7–0.9 represents reasonable/moderate performance; and values > 0.9 indicate high performance (Peterson et al. 2011). The percent variable contribution and jackknife permutation importance functions in Maxent were used to investigate the relative importance of bioclimatic predictors affecting species distributions. The response curves of high contribution values (> 25 %) generated by Maxent were used for investigating the responses of the seven species to environmental variables.

The global model was used to predict the potential distribution of exotic species in Japan for the following reason: if exotic species have not been present for a sufficient length of time to spread to all suitable areas, their distributions are limited to areas smaller than their potential distribution areas. Thus, estimation using the present data in invasive areas may potentially underestimate their distribution areas (Le Maitre et al. 2008).

The invasion stages of the exotic species in Japan were evaluated using an approach that was theoretically developed by Gallien et al. (2012) and applied analytically by Kumar et al. (2015). According to this framework, a species would be at quasi-equilibrium if the regional and global models predict higher probabilities (> 0.5) for the species’ presence (i.e. stabilising populations). In contrast, if both models predict lower probabilities (< 0.5) for species presence, these locations may represent population sinks (i.e. sink populations). If species presences cover the global niches (probability > 0.5) but not regional niches (probability < 0.5), this finding suggests colonisation from different sources, including already invaded areas in the regional invaded range (i.e. colonising populations). In contrast, if species presences cover the regional niches (probability > 0.5), but not global niches (probability < 0.5), these populations may be adapting to new environmental conditions (i.e. local adapted populations).

Seven exotic species were reported from 1,516 locations in Japan and 4,135 and 1,953 location data were found for Europe and North America, respectively (Table 1). Armadillidium vulgare, H. danicus, P. scaber and P. pruinosus were reported from more than 100 locations in Japan and A. vulgare had the widest distribution area, from N43°10' to N24°45'. In contrast, A. nasatum, P. dilatatus and P. laevis were found in less than 32 locations (Figs 1–7, Table 1). Most species were rarely reported in the Hokkaido, but P. scaber was found in many locations in the Hokkaido (Fig. 6). Porcellionides pruinosus was reported at the southernmost location of Japan included in the dataset, Iriomotejima Island (N24°16') (Fig. 7). Five of the seven exotic species were first reported in 1927 and the two other species (A. nasatum and P. dilatatus), with narrow distribution areas in Japan, were reported in 1943 (Table 1).

Figure 1.

Potential distribution maps, invasive stage and occurrence locations for A. nasatum. Potential distribution maps predicted by the regional model (A) and the global model (B). Invasive stage of A. nasatum populations based on global and regional model predictions (C) and the locations on the map (D); red: stabilising population, yellow: sink population. Four large islands and three regions were described in A and B, respectively.

Figure 2.

Potential distribution maps, invasive stage and occurrence locations for A. vulgare. Potential distribution maps predicted by the regional model (A) and the global model (B). Invasive stage of A. vulgare populations based on global and regional model predictions (C) and the locations on the map (D); red: stabilising population, blue: colonising population, yellow: sink population.

Figure 3.

Potential distribution maps, invasive stage and occurrence locations for H. danicus. Potential distribution maps predicted by the regional model (A) and the global model (B). Invasive stage of H. danicus populations based on global and regional model predictions (C) and the locations on the map (D); red: stabilising population, blue: colonising population, green: adapted population, yellow: sink population.

Figure 4.

Potential distribution maps, invasive stage and occurrence locations for P. dilatatus. Potential distribution maps predicted by the regional model (A) and the global model (B). Invasive stage of P. dilatatus populations based on global and regional model predictions (C) and the locations on the map (D); blue: colonising population, green: adapted population.

Figure 5.

Potential distribution maps, invasive stage and occurrence locations for P. laevis. Potential distribution maps predicted by the regional model (A) and the global model (B). Invasive stage of P. laevis populations based on global and regional model predictions (C) and the locations on the map (D); red: stabilising population, blue: colonising population, yellow: sink population.

Figure 6.

Potential distribution maps, invasive stage and occurrence locations for P. scaber. Potential distribution maps predicted by the regional model (A) and the global model (B). Invasive stage of P. scaber populations based on global and regional model predictions (C) and the locations on the map (D); red: stabilising population, blue: colonising population, yellow: sink population.

Figure 7.

Potential distribution maps, invasive stage and occurrence locations for P. pruinosus. Potential distribution maps predicted by the regional model (A) and the global model (B). Invasive stage of P. pruinosus populations based on global and regional model predictions (C) and the locations on the map (D); red: stabilising population, blue: colonising population, yellow: sink population.

Training AUC value for each species varied between the global model (0.681–0.910) and the regional model (0.712–0.857); the values for three (A. nasatum, A. vulgare and P. dilatatus) of the seven species were higher in the global models than in the regional models (Table 2, Suppl. material 5). In the global model, annual mean temperature and mean diurnal temperature range (annual mean of all the month’s maximum - minimum diurnal temperatures) became important predictors (> 25 %) for five species and temperature seasonality also became an important predictor for P. dilatatus (37.6 % contribution). However, the precipitation-related variable, precipitation seasonality was only an important predictor for A. nasatum (25.7 % contribution), P. dilatatus (25.3 % contribution) and P. laevis (28.0 % importance). As well as in the global model, mean diurnal range also became an important predictor (> 25%) for five species in the regional model. However, the annual mean temperature exceeded 25% importance only for A. nasatum (32.2 % importance) and P. pruinosus (39.5 % contribution). In addition to the two climatic variables, maximum temperatures of the warmest quarter (A. vulgare, 25.3 % importance; P. laevis, > 50 % contribution and importance) and the wettest quarter (P. dilatatus, 27.2 % contribution) became important predictors. Within precipitation related-variables, precipitation of the wettest month was only an important predictor for A. nasatum (> 65 % contribution and importance).

In the global model, relationships between high-contribution (> 25 %) variables based on the percent contribution values and the probabilities of seven species are shown in Figs 8 and 9. The relationship between the probability of occurrence and mean diurnal range showed similar trends amongst five species, A. vulgare, P. dilatatus, P. laevis, P. scaber and P. pruinosus; their probabilities of occurrences were high in mean diurnal ranges of less than 10°C. However, there were different trends in the relationships between occurrence probability and annual mean temperature between species; the probability of occurrence of A. nasatum, H. danicus and P. scaber preferred high temperatures (> 15 °C) and P. pruinosus tended to be distributed in areas with more moderate temperatures (15–25 °C). Armadillidium nasatum was associated with low precipitation seasonality. Porcellio dilatatus tended to be distributed with lower temperature seasonality, although this relationship was weak.

Figure 8.

Relationships between high-contribution variables (more than 25 % contribution) and the occurrence probabilities. The contribution variables are based on the percent contribution values in the global model. First row: A. nasatum; second row: A. vulgare; third row: H. danicus; fourth row: P. dilatatus. Left: the highest contribution variable; right: the second highest contribution variable.

Figure 9.

Relationships between high-contribution variables (more than 25 % contribution) and the occurrence probabilities. The contribution variables are based on the percent contribution values in the global model. First row: P. laevis; second row: P. scaber; third row: P. pruinisus. Left: the highest contribution variable; right: the second highest contribution variable.

The potential distributions of the seven species in Japan based on the regional and global models are shown in Figs 1–7. The global models for six species represented wider ranges than the regional ranges, excluding P. dilatatus with a wider potential distribution area in the regional model. The potential distribution maps of the global model for A. nasatum, A. vulgare, H. danicus, P. scaber, and P. pruinosus were almost identical (Figs 1–3, 6 and 7), which indicates that suitable regions would occur in western Honshu, Shikoku and Kyushu. However, there were small differences in the potential distribution ranges amongst the species; those of A. vulgare, P. scaber and P. pruinosus were further inland compared with those of A. nasatum and H. danicus and the potential distribution area of P. scaber was covered coastal areas of Hokkaido. The potential distribution area of P. laevis was limited to more south-western regions than those five species (Fig. 5). Porcellio dilatatus was potentially distributed in most areas of Kyushu and Shikoku and in scattered areas of central Honshu (Fig. 4). Based on comparison with the global model, the potential distribution maps of the regional model showed large interspecific variations. The potential distribution areas of A. vulgare, H. danicus, P. scaber and P. pruinosus were limited to narrow areas within coastal regions (Figs 2, 3, 6 and 7). For A. nasatum, the potential distribution area was large and included a large portion of eastern Honshu and coastal areas of western Honshu (Fig. 1). Only P. dilatatus showed distribution in larger areas in the regional model than in the global model (Fig. 4); in addition to most of Kyushu and Shikoku, as in the global model, the potential distribution area also extended to coastal areas of Honshu and southern Hokkaido. The regional model predicted distribution of P. laevis in the Setouchi, Kanto and Okinawa regions (Fig. 5).

Plotting the probabilities of the regional model against those of the global model showed three patterns (Figs 1–7). The first was that the global probabilities were more than 0.5 regardless of the regional probabilities, which indicated that these species were composed of stabilising and colonising populations, with several locations considered to represent sink populations (less than 0.5 probabilities in both models) also included. This pattern was found for A. vulgare, H. danicus, P. laevis, P. scaber and P. pruinosus (Figs 2, 3 and 5–7). The second pattern was found for A. nasatum and was composed of stabilising populations with one sink population (Fig. 1). The third pattern was shown for P. dilatatus; the probabilities of the regional model were high, whereas those predicted by the global model were low (Fig. 4). These populations may have adapted to new environments (local adapted populations). The four types of populations, i.e. stabilising, colonising, local adapted and sink populations, did not show distinctive geographic patterns, but sink populations tended to be found in Hokkaido and inland on Honshu.

To estimate the factors limiting distribution areas at a macro scale, ecological niche models (ENMs) and species distribution models (SDMs) have been useful (Elith and Leathwick 2009, Franklin and Miller 2010, Peterson et al. 2011, Franklin 2013). However, few studies have applied these methods to the analysis of soil animals, excluding ants (Fitzpatrick et al. 2007, Kumar et al. 2015). The present study was a first attempt to apply ENMs to evaluate the important factors determining the distribution areas and invasion stages of exotic terrestrial isopods. Needless to say, however, the occurrence data used in the study could not cover all distribution areas of the species. Thus, it is possible that the ENMs constructed will be altered when a more accurate and larger database has been developed.

The global model indicated that temperature-related variables are more important than precipitation-related variables in limiting the distributions of terrestrial isopods. The distributions of isopods at a macro scale are limited by natural factors, especially temperature and moisture (Harding and Sutton 1985, Hopkin 1991). For example, a cartographic analysis of terrestrial isopods in the former USSR indicated that the length of the period with temperatures above 10 °C plays an important role in limiting isopod distributions (Kuznetsova and Gongalsky 2012). Lilleskov et al. (2008) also found that the species richness of exotic terrestrial isopods in North America peaked at mid-latitudes (N40°). These facts supported the inference that temperature may be an important factor limiting the distribution of exotic terrestrial isopods in Japan. This study also found that, not only annual mean temperatures but also mean diurnal temperature ranges, played important roles in delimiting the distributions of exotic species. Armadillidium vulgare, P. dilatatus, P. laevis, P. scaber and P. pruinosus were potentially distributed in areas where diurnal temperature fluctuated less. However, there were differences in suitable annual mean temperatures between species. Armadillidium nasatum, H. danicus and P. scaber preferred warm regions, whereas the suitable temperature range of P. pruinosus showed a hump-shaped pattern. However, the annual mean temperatures of high probability (probability greater 0.5) for these species ranged from 10 °C to 25 °C, which covered most areas of Japan, excluding mountainous regions. This result suggested that the annual mean temperature was not an important factor limiting the distributions of these species in Japan and that climate warming may cause latitudinal shifts in distribution ranges for terrestrial isopods (Hickling et al. 2006, Mason et al. 2015) and there is a possibility that their ecological functions may change. However, there is no room for doubt that local-scale environments, e.g. microbial activity and vegetation, play important roles in limiting distributions of terrestrial isopods (Zimmer and Topp 1999, 2000). Thus, an integrated analysis of local- and macro-scale features will be required to understand the factors limiting the species distribution (Pateman et al 2016).

For six of the seven exotic isopod species, the global model predicted wider potential distribution areas in Japan than did the regional model. These results implied that these species have not yet filled all suitable environments in Japan and that there is a risk that their distribution areas will spread further in the future. The results also indicated that combined data from native and introduced regions were useful for estimating the potential distributions of these exotic species in the invaded region (Le Maitre et al. 2008, Elith 2015, Ray et al. 2016). The present study evaluated the invasive stages of these species in Japan by applying the theoretical framework of Gallien et al. (2012), which has been successfully used to evaluate the invasive stages of some exotic species (Kumar et al. 2015, Zhu et al. 2017). The evaluation revealed that many populations of five species (A. vulgare, H. danicus, P. laevis, P. scaber and P. pruinosus) in Japan were stabilising and colonising. A stabilising population was defined as a population predicted by both models as high probability (Gallien et al. 2012, Kumar et al. 2015) and indicates that these populations may be maintained continuously. A colonising population was a population distributed in an area where the regional model predicted low probability, although these areas were predicted as high probability by the global model (Gallien et al. 2012, Kumar et al. 2015). The occurrence of many colonising populations suggests potential for more populations to be established in the area, which will be changed to high probability in the regional model as well. Thus, evaluating the invasive stage also enables prediction of the future spread of distribution areas of these five species. In addition, the species’ populations included several sink populations distributed in areas characterised as low probability in both models (Gallien et al. 2012, Kumar et al. 2015). These populations were found inland on Honshu and Hokkaido, which may have been too cool to maintain these populations.

Armadillidium nasatum and P. dilatatus were reported from relatively few locations in Japan (20 and 7 locations, respectively). Examination of the literature revealed that A. nasatum and P. dilatatus were first reported 16 years after the earliest reports of the other species. It is speculated that shortage of elapsed time since the introduction of A. nasatum was one of reasons for its distribution in fewer locations in Japan, because the species’ suitable area covered a large area of Japan and because this species had been introduced to suitable areas; this species was composed of stabilising populations. Additionally new locations were recently found (Yamasaki and Yoshida 2017, T. Suzuki, personal comments); thus, it is possible that this species will spread its distribution areas in the future, in addition to the five species mentioned above. However, for the narrow distribution range of P. dilatatus, the ecological niches and introduced areas of this species were more important for explaining its occurrences in fewer locations rather than the time elapsed. The relationship between the regional and global models for P. dilatatus presented different traits from those of the other species. The regional model for P. dilatatus predicted the south-western region and coastal area of Honshu as suitable areas, but the global model indicated Shikoku and Kyushu as suitable areas. These results suggested that this species was introduced in unsuitable environments and adapted to these new environments in Japan, which may be a reason why this species was found in fewer locations. In addition, the global model indicated that this species could spread if it migrates to south-western areas (i.e. Shikoku and Kyushu). However, it is possible that the small number of occurrence locations in Japan was insufficient to precisely estimate the potential distribution map (Hernandez et al. 2006, Wisz et al. 2008). Therefore, further research is required to evaluate the adaptations of these taxa to new environments in Japan.