The study region includes most parts of the European mainland and the larger islands (excluding only European Russia, Ukraine, Belarus and Cyprus) (Fig. 1). To provide an up-to-date overview of the current distribution of the coypu, data from several sources such as publications, national administrative authorities and scientists were compiled between July 2017 and December 2018 (Suppl. material 1, Table S1). Additionally, we downloaded occurrence data from the Global Biodiversity Information Facility (GBIF.org). The occurrence of persisting populations in the Boreal, Arctic and northern Alpine biogeographic regions (European Environment Agency 2002) is not mentioned in literature (Carter and Leonard 2002, DAISIE 2009, Tsiamis et al. 2017), thus, no further effort has been made in sending out data requests for those regions. The same applies to Portugal, from where occurrences are not known either (Carter and Leonard 2002, Tsiamis et al. 2017). We did not include data from the coypu’s native range, because there are only few spatially explicit records available, as already pointed out by Jarnevich et al. (2017).

Figure 1.

Coypu occurrence records from 1980 to 2018 in Europe. The decade-wise accumulation of records is depicted on the left side, with records of the respective decade in black and records of previous decades shown in grey. On the summary map (right side), records in Great Britain are displayed in red, as the coypu is officially eradicated (see Gosling and Baker 1989). Note, that for optimised illustration purposes, 10-km buffered centroids of occurrence records are shown.

The raw occurrence data were prepared and quality-checked prior to analyses. If a source described the occurrence of coypu over several decades, the record was split into one record per decade. Only records that contained geographic information and approximate sampling date were considered in this study. Records that were lacking coordinates, but contained an unambiguous locality description, were georeferenced, either using the point-radius method (estimating coordinates and an uncertainty radius according to the precision of the locality description) or the shape method (assigning a geographic shape that represents the uncertainty) (Wieczorek et al. 2004); i.e. for France and Germany, locality descriptions at municipality levels were linked to the according feature of the GADM 3.4 shapefile (database of Global Administrative Areas, https://gadm.org/). Point-radius georeferencing was conducted using the GeoLocate web tool (Rios and Bart 2010) and web map services. Coordinate uncertainty estimates were used to buffer the records. Note that for records linked with the municipality area, no buffers were introduced; thus, the uncertainty information was linked with the shape and size of the municipalities. As the coypu is a mobile species, a circular buffer of 1 km radius was applied to records with very low uncertainty estimates to sufficiently cover potential home ranges and account for cases were records fall at the borders of grid cells. After a literature research on coypu home ranges (Doncaster and Micol 1989, Denena et al. 2003, Nolfo-Clements 2009, Scheide 2013), we have orientated ourselves to the upper end of documented values. Scheide (2013) reported that, along waterways, territory length usually varies between around 150 m to about 1 km, although home range can increase up to several kilometres in radius if resources are scarce.

For further analysis, records that were missing essential information (i.e. no georeference), putative duplicates (i.e. records with same year and coordinates or locality description), as well as records exceeding an uncertainty radius of 10 km in areas where more accurate records were available, were discarded. This resulted in a final dataset consisting of 24,232 coypu records between 1980 and 2018 across 28 European countries containing year, uncertainty estimate and coordinates (Fig. 1).

We transformed those presence records (inclusive uncertainty buffer) to a grid of 5 × 5 km resolution (temporal resolution: one year) to reduce the effect of pseudo-replication (i.e. artificial inflation of the sample size due to intensively-sampled regions or non-detected duplicates in different datasets). Grid cells that showed only marginal overlap with buffered presences (< 2.5% of the grid cell area) were not defined as presences, to avoid area inflation. After discarding those, 25,534 grid cells (about 12.6% of all grid cells throughout the study area) were defined as presence grid cells and were used for further analysis (see workflow scheme, Fig. 2). Note, that there are more presence grid cells as presence records per se, as we used the buffered presences for the grid cell transformation to consider the spatial uncertainty of records.

Figure 2.

Workflow used for species distribution modelling. Presence data was transformed to a raster grid (5 × 5 km), spatiotemporally analysed, accordingly grouped into five levels of persistence and then combined with environmental data and generated absences (using four different sampling strategies) to fit logistic regression models. The absence sampling strategy that derived the best evaluation measures was used for predicting suitability under current (1979–2013) and future climatic conditions after going through a model selection and averaging procedure.

Dullinger et al. (2009) found that environmental niche models, based on distribution data, produce more accurate predictions when analyses are restricted to persistent populations. To keep information about long term occurrence and, therefore, probable establishment, we conducted a spatiotemporal analysis and allocated the presence within grid cells to different persistence categories. Presence records per grid cell were analysed by counting the number of years containing coypu records and the time-span covered. Grid cells with multiple years of recorded coypu presence were classified with regard to the time-span covered by the records and accordingly grouped into different persistence levels of at least one, at least two or at least three generation lengths (hereafter GL; which is 5 years, according to Ojeda et al. (2017)). Temporal coverage of at least one or two generation lengths can be interpreted as early signs of establishment. For cells that repeatedly showed coypu presence, covering a time-span of at least 15 years (thus three full generation lengths of coypu), we assumed the occurrence of established populations within the time from 1980 to 2018. If temporal discontinuities of more than 10 years between successive records were present in a grid cell, we divided the counts into subgroups, which were analysed separately. In this case, the highest derived persistence category was assigned to the according grid cell.

On this basis, we created sub-datasets by stepwise exclusion of the lowest level of persistence (those sub-datasets are hereafter called persistence levels and abbreviated as indicated by the bold words: all grid cells = 25,534, multiple records per grid cell = 15,078, multiple records per grid cell that cover at least one GL = 11,619, multiple records per grid cell that cover at least two GL = 5,145, established populations only = 2,505).

Environmental predictors included bioclimatic variables from the CHELSA database (Karger et al. 2017a), CORINE Land Cover data for different land cover classes (i.e. agricultural areas, wetlands and water bodies) (European Environment Agency (2012), hilliness (i.e. the standard deviation of the average sea level) derived from a digital elevation model (European Environment Agency 2000), human population density (EUROSTAT 2011) and data on rivers and lakes (European Environment Agency 2012b) (Table 1). Human population density was log-transformed to reduce the effect of outliers due to large cities. All variables were calculated for a grid of 5 × 5 km and standardised for better comparison (i.e. scaled to a mean of 0 and a standard deviation of 1). As multiple bioclimatic variables were highly correlated (|Pearson’s r| > 0.7), we limited them to five variables that represent temperature (bio06; ‘Minimum Temperature of Coldest Month’ and bio10; ‘Mean Temperature of Warmest Quarter’), temperature fluctuations (bio2; ‘Mean Diurnal Range’), precipitation (bio17; ‘Precipitation of Driest Quarter’) and precipitation fluctuations (bio15; ‘Precipitation seasonality’). After this procedure, all pairwise correlations of our final predictor set were below the threshold of |Pearson’s r| > 0.7 and all variables showed VIF (Variance Inflation Factor) values smaller than 10, indicating that collinearity is a minor issue within the predictor data (R package usdm, Naimi et al. 2014) (for correlations amongst variables, see Suppl. material 1, Fig. S1).

Additionally, we used qualitative information on biogeographic region (European Environment Agency 2016) and country (GADM 2018) for summary statistics.

Two IPCC (Intergovernmental Panel on Climate Change) climate change scenarios were selected to model coypu response to a changing climate by the mid-21^st century (2041–2060). One represents medium (RCP 4.5; Representative Concentration Pathway) and one represents severe climate change (RCP 8.5) by depicting the different approximate radiative forcing in comparison to the pre-industrial state (i.e. + 4.5 and + 8.5 W/m²) (Moss et al. 2010). Further, it is known that climate predictions are sensitive to different climate model frameworks (GCMs; General Circulation Models) (Porfirio et al. 2014). Thus, we downloaded data for two different GCMs, representative for different model families (Sanderson et al. 2015), from the CHELSA website (Karger et al. 2017a). The GCMs were chosen considering model independence and performance: Had-GEM2-A0 and CESM1-BGC.

To model the range of potential distribution under current climate and under climate change, we used logistic regression, a generalised linear modelling (GLM) technique that is widely used for predicting species distributions (Elith and Leathwick 2009, Franklin 2010). Here, a linear model is related to the binary response variable via a logistic link function (McCullagh and Nelder 1989). We chose this technique due to the proposed transferability in time and space and lower risk of overfitting compared to other methods, such as classification trees (Marmion et al. 2009, Franklin 2010). Since logistic regression requires presence and absence data, we generated absences using different approaches to compare their model performance and predictions in order to select the most appropriate absence design. The absences were drawn either as background or pseudo-absences, i.e. from all grid cells or only from non-presence grid cells (see Phillips et al. 2009); this was done for the whole study area extent or within a buffer of 150 km around presence grids, representing an assumption of coypu dispersal (Aliev 1968, Gosling and Baker 1989, Hong et al. 2015) . As the data collected consists of a variety of sources, contains opportunistic records and is not following a standardised collecting scheme across our study area, we used a target-group approach to account for biased survey effort (Phillips et al. 2009, Stokland et al. 2011). Here, we filtered occurrence data of non-marine, non-volant, small- to medium-sized European mammals from GBIF (https://www.gbif.org/) for the same spatial and temporal extent and downloaded the resulting records, assuming these will exhibit a similar spatial sampling bias. The records were used to create a probability density surface (Suppl. material 1, Fig. S2a) that served as the basis for generating absences.

We combined each of the four different absence sampling strategies (background full, background restricted, pseudo-absence full, pseudo-absence restricted; see Suppl. material 1, Fig. S2b) with all five levels of persistence, resulting in a total of 20 data designs. The number of absence grids cells was set to at least 10,000 or equal the number of presence grid cells. Barbet-Massin et al. (2012) argue that this amount of generated absences adequately depicts the model quality without the need to account for variability in generated absences, for example, through generating replicates. At the same time, predictive accuracy of GLMs does not significantly increase with prevalence, once the number of presences reached at least one tenth of the number of absences (Barbet-Massin et al. 2012). In their study, this held true for weighted and unweighted schemes, therefore we did not apply any weightings.

As pseudo-absences which were generated across the whole study area extent (“pseudo-absence full”) consistently derived the best evaluation values (Suppl. material 1, Fig. S2c) and spatial predictions were basically identical between different absence sampling strategies, we chose this sampling strategy for further analysis.

We assessed the goodness of fit for the full models (including linear and quadratic terms for all variables) of all datasets and of 25 random subsets per persistence level. Those were created by drawing random presence subsamples that equal the size of the according persistence level to check for sampling size effects that might occur. For the model evaluation, we compared a set of commonly-used measures (Allouche et al. 2006, Elith and Leathwick 2009, Liu et al. 2011), based on a fivefold cross-validation (split ratio train:test equals 80:20): 1) AUC, which is the sum of the area under the receiver operating curve, a graph that displays false positive vs. true positive rate (Franklin 2010); 2) adjusted D², which is the proportion of explained deviance, taking into account the number of model parameters and observations and thus allows comparison amongst different models (R package modEvA; Barbosa et al. 2014); and 3) the threshold-dependent true skill statistic (TSS) which needs a binary result. TSS corresponds to the sum of sensitivity (i.e. the proportion of correctly-predicted presences) and specificity (proportion of correctly-predicted absences) minus one and was shown to be independent of prevalence (Allouche et al. 2006). As we aimed to identify the regions sensitive to coypu invasion and regarded false negatives as costlier than false positives, the threshold for TSS computation was chosen following the recommendation of Jiménez-Valverde et al. (2011) to avoid omission error by maximising sensitivity whilst keeping a reasonable specificity and therefore TSS value. Therefore, we set the sensitivity to a fixed value of 0.95 (i.e. the threshold used for binary classification will lead to 95% of presences predicted correctly).

Variable importance was measured for each predictor by evaluating the mean drop in explained deviance caused by removal of the respective predictor. Finally, true positives and omission errors were mapped to reveal sensitivity issues and spatial patterns in model performance (Suppl. material 1, Fig. S3).

The full models’ quality notably increased with increasing level of persistence and this effect could be clearly distinguished from sample size effects, when comparing the evaluation measures with those of the random subsets (Fig. 3). We hence used the persistence level classification for prediction. To select the best models for each persistence class, we compared the corrected Akaike’s Information Criterion (AICc) between the full model (all predictors) and possible sub-models. Then we averaged the top models with ΔAICc < 4 (R package “MuMIn”, Barton 2019), finally resulting in five averaged models (one per persistence level) that were fitted with the whole data of the according persistence level.

Figure 3.

Comparison of model performance for the persistence levels (with pseudo-absences for the full extent) versus random subsets. The boxplots show the results of the fivefold-cross validation for the persistence levels (blue) and the according random subsets (orange).

We assumed grid cells that contain long-term occurrences to be more informative than those where coypu occurrence was only registered once. Both approaches (using all data vs. subsets) in its extremes may incorporate biases (i.e. all presence data will more likely include non-persistent occurrences and false identifications, whereas grid cells that show long-term occurrence might underestimate the area of a still-spreading alien species and comprise historical effects of propagule pressure due to regional differences in fur farming intensity, as well as effects of uneven data availability across regions). To balance those possible biases and reduce uncertainty, we combined the resulting predictions of probability of occurrence for all persistence levels and created a consensus prediction by simply calculating the mean probability of occurrence per grid cell, depicting the overall agreement of the averaged models. Marmion et al. (2009) showed a significant increase in accuracy and robustness for consensus predictions that used averaging methods. For binary maps, we used the same cut-off as for the computation of the TSS (sensitivity fixed to 95%) to separate suitable and unsuitable grid cells (Jiménez-Valverde et al. 2011).The proportion of suitable grid cells was calculated for all countries and compared to the proportion of grid cells that contained presence records.

Further, consensus forecasts for all climate change scenarios were computed. The change in probability of occurrence was assessed by comparing the number of cells that were classified as suitable under current climate and under climate change scenarios and by subtracting the probabilities of occurrence of current from future predictions. To obtain the agreement between binary models, we used the sum of predicted presences per cell across all averaged models, with a high value meaning high agreement.

Finally, we used the resulting predictions to define priority regions for surveillance and management by creating a risk map. Grid cells that 1) show high probability of occurrence in the consensus prediction and 2) areas adjacent to already known recent occurrences, are deemed to be particularly susceptible to invasion by coypus, due to short colonisation distances. Thus, to incorporate dispersal constraints and to account for proximity to known occurrence, a weighting matrix was computed, by summing up weighted inverse Euclidean distance classes per decade for each cell (Suppl. material 1, Fig. S4).We combined this matrix with the consensus map under current climate showing the mean probability of occurrence. Here, values below the probability threshold for binary map computation were set to zero. We expected both suitability and proximity to be of equal importance for the invasion process. If a cell were considered too distant or unsuitable, no risk of invasion was assumed. For the United Kingdom, we considered only Northern Ireland for the risk map and excluded Great Britain, as the coypu has been eradicated there (see Gosling & Baker 1989).

Statistical analyses were conducted and maps were produced using ArcGIS 10.5.1 (ESRI 2018), R 3.6.0 (R Core Team 2019) within the GUI RStudio 1.1.463 (RStudio Team 2018). The following R packages were used: dismo (Hijmans et al. 2017), dplyr (Wickham et al. 2017), modEva (Barbosa et al. 2016), MuMIn (Barton 2019), PresenceAbsence (Freeman and Moisen 2008), raster (Hijmans 2017), spdep (Bivand et al. 2013, Bivand and Piras 2015), rgbif (Chamberlain et al. 2017), usdm (Naimi et al. 2014).

In total, 24,232 coypu presence records (corresponding to 25,534 grid cells at 5 × 5 km) were collected across 28 countries. The spatiotemporal analysis of presence grid cells shows centres of documented long-term occurrence in Czech Republic, France, Germany, Italy and the Netherlands. Of all presence grid cells, 45.5% (corresponding to 20 countries) show at least early signs of establishment (i.e. had multiple records that covered one generation length as a minimum; of those 20.1% have been covered by at least two generation lengths and 9.8% of the grid cells (corresponding to 10 countries) show spatially-explicit evidence for long-term persistence (i.e. established populations) with coypus being present over a period of at least 15 years (Fig. 4). Note, that these periods of occurrence do not necessarily imply that populations are still present in a given area, but are indicative of the general suitability of the area within the last decades. For example, the successful eradication programme in Great Britain led to the local extinction of the coypu (Gosling and Baker 1989), despite the suitability of the environment that allowed persistence over several generation lengths.

Figure 4.

Persistence levels of presence grid cells as derived from the spatiotemporal analysis. Each grid cell that intersects at least one record of coypu presence between 1980 and 2018 is coloured according to the maximum derived persistence level: 1) single record, 2) multiple records, 3) one generation length (multiple records covering at least 5 years), 4) two generation lengths (multiple records covering at least 10 years), 5) established (multiple records covering at least 15 years). One generation length is assumed to be 5 years, following Ojeda et al. (2016).

Model quality increased with higher levels of persistence, with mean AUC values ranging from 0.90 (all) to 0.96 (established) indicating excellent discrimination ability across all averaged models and TSS values ranging from 0.61 to 0.79 which can be interpreted as good to excellent agreement between training and test data (Eskildsen et al. 2013). Adjusted explained deviance was between 43.0% and 57.5% (Table 2).

Between two to four top models were averaged for the persistence levels, with the full model always being included. Only land cover variables were excluded (‘shores’, ‘water-bodies’) and none of them was excluded across all persistence levels (for model weightings and ΔAICc, see Suppl. material 1, Table S2). The analysis of predictor variable importance showed that the mean temperature of the warmest quarter (bio10; mean drop in D² ± SD: 8.5 ± 0.77), mean diurnal temperature range (bio2; 7.8 ± 3.51) and the minimum temperature of the coldest month (bio6; 3.2 ± 0.22) were the most important of the analysed predictors (Fig. 5). In addition, precipitation seasonality (bio15; 2.1 ± 0.41) played a relatively important role, whereas the other variables had markedly lower values. Still, the mean precipitation of the driest quarter (bio17), hilliness and the distance to settlements are of superior importance in comparison to the land cover variables. Coypu presence was more likely at medium diurnal temperature ranges and when the mean temperature of the warmest quarter, as well as the minimum temperature of the coldest month, was medium to high. The probability of occurrence decreased with increasing precipitation seasonality, distance to settlements and hilliness and increased with increasing human population density, number of wetlands and higher precipitation during drier months (bio17).

Figure 5.

Importance of predictor variables for the averaged final models. The importance is measured as the mean drop in explained deviance (D²) upon removal of the respective predictor. For descriptions of the predictor variables, see Table 1.

Plotting of the omission errors of the binary consensus prediction revealed that those mostly occurred at the range margins of the currently known European distribution, especially towards Southern Europe and mountainous areas (Suppl. material 1, Fig. S3b). A good distinction between the presence and generated pseudo-absence data could be achieved (median of predicted probability for pseudo-absences and presences: 0.05 vs. 0.77, Figure S3a).

The consensus map for current climatic conditions shows that, currently, large parts of Europe have a high probability of coypu occurrence (Fig. 6a). Applying a threshold that gives 95% sensitivity (correctly predicted presence grid cells) results in 42.9% of the study area being rated as potentially suitable (Fig. 6b). Only 27.4% of those cells already comprise documented coypu occurrences, while the remaining 72.6% being, to our current knowledge, yet unoccupied. These potentially suitable areas cover most of Central Europe and parts of the following biogeographical regions; Atlantic (67.8%), Black Sea (92.7%), Continental (79.3%), Mediterranean (24.0%), Pannonian (93.5%) and Steppic (38.2%) regions. Only minor parts of the Alpine and Boreal regions contain predicted suitable grid cells (8.8% and 4.2%) and the Arctic biogeographic region is considered unsuitable.

Figure 6.

Consensus predictions of the probability of occurrence across the study area under current climatic conditions (years 1979–2013). a Mean probability of occurrence across the final averaged models of all persistence levels. b Binary classification of suitable and unsuitable grid cells after applying a threshold corresponding to 0.95 sensitivity (= 0.16).

All four climate change scenarios show substantial shifts in predicted habitat suitability until the mid-21^st century (2041–2060) (Fig. 7). While the total amount of suitable area is predicted to decrease between 2–8% in comparison to current climatic conditions (38.1% (HadGEM1 A0 RCP 4.5), 34.7% (HadGEM1-A0 RCP 8.5), 40.9% (CESM1-BGC RCP 4.5) and 39.8% (CESM1-BGC RCP 8.5), particularly northern and Atlantic regions with Ireland and the United Kingdom will experience an increase in suitability. In contrast, all models predict decreasing suitability along the southern range.

Figure 7.

Future predictions. Agreement between averaged models for projected probability of occurrence in the mid-21st century under two climate change scenarios (medium climate change: RCP 4.5; severe climate change: RCP 8.5) combined with two different global circulation models (HadGEM1-A0, CESM1-BGC, displayed as number of models predicting presence (left side) and net change in occurrence probability compared to the current climatic situation (right side).

The risk map (Fig. 8) is an indication of invasion risk, considering not only the potential suitability, but also the current distribution of the coypu and thus the likelihood of dispersal and colonisation of new grid cells. In Figure 9, the percentage of suitable and occupied area per country, as well as the establishment status of the coypu in the respective country, is shown. None of the already affected European countries reached saturation by now and, additionally, a number of not-yet invaded countries contain a considerable amount of suitable area.

Figure 8.

Risk map, highlighting regions potentially prone to invasion, i.e. with high probability of occurrence under current climate and adjacent to known recent occurrences of coypu. Presence grid cells are shown in black. Great Britain was excluded as the coypu is officially eradicated (Gosling & Baker, 1989).

This study confirms and substantially expands the overview of Tsiamis et al. (2017), who report coypu occurrence for 18 EU countries and establishment for 12 of those. Our presence dataset covers 28 European countries, of which 10 countries showed spatially-explicit evidence for long-term establishment of the species and another 10 countries showed at least early signs of establishment within the regarded time period (Fig. 9; note that our study also deals with non-EU member states). Under the current climate, a considerable number of countries has high proportions of suitable areas, for example, Belgium, Germany, Luxembourg, Netherlands and Hungary with > 90% and Bulgaria, Croatia, Czech Republic, Denmark, France, Kosovo, Poland, Italy, Serbia, Slovenia, Slovakia and the United Kingdom with more than half of their country area being predicted as suitable for coypu (Fig. 9). The comparison of presence grid cells with not-yet invaded but suitable ones, shows that further substantial range expansions can be expected. Moreover, several countries that do not have documented coypu occurrences yet, contain potentially suitable areas. Overall, 42.9% of the study area is considered suitable under current climate (1979–2013) (Fig. 6), of which less than a third already contained occurrences from 1980 to 2018.

Figure 9.

The country-wise percentage of suitable grid cells under current climate (brown bars) and grid cells containing presences (red bars). Countries are marked according to the maximum persistence level. Countries with no suitable areas and occurrences are not shown, as well as microstates not covering a whole grid cell. Alphabetically ordered country abbreviations with corresponding percentages of suitable grid cells and occupied cells in parentheses: AL : Albania (4.7/0.2), AT : Austria (30.0/6.3), BA : Bosnia and Herzegovina (48.2/-), BE : Belgium (99.9/15.6), BG : Bulgaria (63.3/2.8), CH : Switzerland (43.2/3.9), CZ : Czech Republic (50.5/26.0), DE : Germany (95.7/61.0), DK : Denmark (70.1/7.7), ES : Spain (17.3/0.6), FI : Finland (-/0.1), FR : France (89.7/55.2), GR : Greece (5.5/2.8), HR : Croatia (89.6/0.8), HU : Hungary (93.6/1.8), IE : Ireland (39.5/2.5), IT : Italy (55.4/10.0), LI : Liechtenstein (33.3/-), LT : Lithuania (29.3/-), LU : Luxembourg (100/19.8), LV : Latvia (4.5/-), ME : Montenegro (14.7/1.1), MK : Macedonia (16.5/2.2), NL : Netherlands (98.3/35.5), NO : Norway (0.8/-), PL : Poland (79.8/0.7), PT : Portugal (12.6/-), RO : Romania (38.6/0.2), RS : Serbia (63.2/0.5), SE : Sweden (6.0/<0.1), SI : Slovenia (74.1/4.1), SK : Slovakia (53.2/6.1), TR : Turkey* (32.4/2.8), UK : United Kingdom (50.3/6.4), XK : Kosovo (64.2/-).*) only the area of the European part of Turkey is considered.

All four climate change scenarios used in this study predicted a slight to moderate decrease of suitable area (from 42.9% under current climate to between 34.7% and 40.9%). This decline is caused by a loss of suitable habitats in the southern parts of Europe, which is not fully compensated for by increasing suitability at higher latitudes (Fig. 7). Thus, our results show that climate change likely will not cause an overall increase of suitable areas for coypu in Europe. This is in line with Bellard et al. (2018) who reviewed modelling studies of climate change effects on alien species distributions and found that climate change will more frequently contribute to a decrease in alien vertebrate species range size. Currently suitable areas closely match warm temperate climates (fully humid or summer dry and with warm summers) after the Köppen-Geiger Climate Classification (Kottek et al. 2006), but not regions with hot summers. Scheide (2013) mentions increasing mortality rates of coypu at high temperatures, which would be in accordance with our finding of decreasing suitability in warming arid regions. In addition, Jarnevich et al. (2017) found support for upper thermal tolerance thresholds of the species. Decreasing suitability on the Iberian peninsula supports the findings of Gallardo and Capdevila (2018) which conducted a risk analysis for Spanish national parks using climate scenarios for 2050 and 2070; they predicted slight to medium decrease of suitability for the coypu in the majority of cases.

The recent occurrence of the coypu in Ireland caused the first Species Alert issued by a European Union Member State under the EU Regulation on Invasive Alien Species. Although those areas were only partly predicted correctly, a considerable area of Ireland is classified as suitable by our predictions. In line with our study, Scheide (2013) and Jarnevich et al. (2017) identified Ireland as having a high similarity to the coypu’s realised niche when predicting its potential distribution on a global scale. Therefore, the recent occurrence reports should be taken with great caution, especially as the overall suitability of Ireland under climate change scenarios is expected to increase (Fig. 7). This said, the predictions of Jarnevich et al. (2017), when modelling the potential distribution of the coypu in the US and worldwide, did also classify large parts of Central and Eastern Europe as unsuitable, which is in contrast with our findings and with the presence of established populations in many areas of those regions.

While our predictions classify most of the Atlantic, Continental, Black Sea and Pannonian biogeographic regions as suitable, this is not the case for the Alpine biogeographic region. Therefore, this study is in agreement with others that have classified the coypu as a typical lowland species and implies that mountain regions act as effective dispersal barriers on a regional scale (Woods et al. 1992, DVWK 1997, Bertolino and Ingegno 2009, Scheide 2013).

Our results highlight the importance of temperature-related climatic variables, such as the mean temperature of the warmest quarter, the mean diurnal temperature range and the minimum temperature of the coldest month as being essential in shaping habitat suitability for coypus (Fig. 5). Under climate change, increasing populations due to decreasing winter mortality seem to be possible and could have economic and environmental consequences in affected areas (Gosling and Baker 1989, Carter and Leonard 2002, Scheide 2013).

Currently, urban coypu populations, fed by humans and profiting from mild urban climate and, in some cases, the thermal pollution of rivers, clearly demonstrate the consequences of high reproduction rates coupled with lowered mortality and enhanced resource availability for population densities (Carter and Leonard 2002, Verbeylen 2002, Walther et al. 2011, Scheide 2013). Whereas Meyer (2005) found local adaptations of the coypu to urban areas, other studies revealed negative effects of settlements on coypu occurrence (Bertolino and Ingegno 2009). In our study, increasing distance to settlements had a negative effect on occurrence probability, whereas the human population density was slightly positively correlated throughout all models. These results indicate that human presence seems to favour coypu occurrence, as it can well adapt to urban waters and can take advantage of additional resources provided for feeding. However, recording bias may contribute to this result as there may be preferential recording in more densely populated regions. Nevertheless, because the attitude toward coypu varies between regions or countries (Carter and Leonard 2002), the association between humans and coypu occurrence may vary spatially and may also change over time. While in some regions the species is hunted (Carter and Leonard 2002, Bertolino and Ingegno 2009) and therefore presence might be more likely in remote areas, in others, regularly fed urban populations occur (e.g. in Germany (Scheide 2013), the Czech Republic (M Anděra personal communication) or in Austria (A Schertler personal observation)). In contrast to former findings, nowadays the coypu is also common in Italian cities, where it is fed by people, highlighting the temporal dynamic of its relationship with humans (S Bertolino personal communication).

Although on a continental scale, climatic aspects are clearly of higher importance (Franklin 2010), some of the predictors associated with land cover consistently showed a significant influence on occurrence probability which, for example, increased with increasing amount of wetlands within a grid cell. The relevance of land cover variables for coypu occurrence was reported by previous studies (Scheide 2013, Farashi and Najafabadi 2015) and should definitely be taken into account in future studies that are conducted on a finer scale. Bertolino and Ingegno (2009) showed that coypu prefers rice paddies as habitats in Northern Italy. Specific agricultural areas do not only enhance food availability, but also potentially provide habitat, for example, trough irrigation ditches. The flexibility to colonise a variety of habitats must be considered when predicting the potential future distribution of the coypu. Although arid environments in Southern Europe are predicted to be unsuitable in the near future (Fig. 7), wetlands and small patches of suitable habitat due to microclimatic factors (e.g. along riparian areas), as well as agroecosystems can still provide suitable habitat for coypus and simultaneously promote conflicts due to feeding damage.

The assumption of an equilibrium between a population and its environment is typically violated during biological invasions, due to ongoing dispersal. Václavík and Meentemeyer (2012) found that SDMs calibrated in early invasion stages tend to be less accurate and under-predict potential ranges of species. The coypu was introduced to Europe more than a century ago and a multitude of release or escape events across regions happened, as it was widely used for fur farming (Carter and Leonard 2002). As the coypu is a conspicuous species that is regularly recorded by a wide range of people (e.g. naturalists, waterway authorities, anglers, farmers and hunters) and given the exhaustive search of records performed, we are convinced that the collated distribution dataset closely reflects the known distribution of the coypu in Europe. Therefore, we assume that our dataset captures a wide range of suitable environmental conditions, with the consensus predictions providing valuable tools to predict the next phase of invasion and areas at high risk (Fig. 8).

The SDMs performed well (Table 2), although some tendency in misclassifications was detected (Suppl. material 1, Fig. S3b). This might be due to missing essential predictors, model mis-specification or influential spatially clustered factors, such as biotic interactions, propagule pressure and dispersal (Elith and Leathwick 2009). Propagule pressure and the number of release events within the last century are difficult to reconstruct given the lack of necessary data, but were presumably differing across Europe, due to the varying economic importance of the fur industry amongst countries (Carter and Leonard 2002). We aimed to account for uneven sampling effort across regions and unsuccessful escape events by using a target-group approach for absence generating and by spatiotemporally analysing presence grid cells. Nevertheless, several escape events at a given site might lead to overestimation of persistence, whereas other regions due to scarcity of spatially explicit records might be under-represented and their suitability hence underestimated (e.g. Southern Europe). Due to these data limitations and because the coypu most likely has not yet colonised all climatically suitable regions in Europe, here, the calculated environmental niche might be conservative. Misclassification increased with increasing hilliness in a region, likely due the chosen spatial resolution, which is too coarse to properly characterise the full variation of environmental conditions in heterogeneous areas (e.g. valleys, which can locally provide suitable habitat in mountainous regions). In addition, a temporal change in sampling effort could lead to earlier decades being under-represented, hence, the impression of the species’ spread could be intensified by more data becoming available recently, for example, through citizen science projects.

The majority of grid cells deemed suitable for coypu under current climate or climate change are not yet colonised. Our results illustrate the urgent need to not only improve management measures in areas with persisting populations, but also find strategies to prevent or reduce further spread as the costs of early intervention are much smaller than control of established populations (Panzacchi et al. 2007). Although the eradication of the coypu in Great Britain was successful, this was achieved as a result of coordinated intensive trapping efforts which were executed by employed professionals over one full decade (Gosling and Baker 1989). Moreover, Great Britain has the advantage of being an island and re-invasion is therefore unlikely. Considering the current situation on the European mainland, with widespread occurrence in Central Europe, it seems highly unrealistic to attempt total eradication of the species in Europe.

Baroch and Hafner (2002) argue that, in the case of low population densities, impacts of coypus in general are rather minor. Given that the coypu is already fairly widely distributed, management that minimises population density and therefore negative economic and environmental impacts should be the aim. As there are several hotspots of coypu occurrence covering more than one country, there is a need for international collaboration to coordinate control measures on a metapopulation scale and prevent compensatory re-invasion from adjoining populations (Oliver et al. 2016). There are well-known cases of migration events of coypus from neighbouring countries, for example, from France to Germany (DVWK 1997, Scheide 2013) and from Belgium and Germany to the Netherlands (Carter and Leonard 2002). Already existing binational control programmes could serve as best-practice examples. Gosling and Baker (1989) suggest concentrating efforts on high-density hotspots to maximise mortality and minimise dispersal to new habitats. Further research regarding coypu dispersal movement and interaction on metapopulation level, also taking into account population genetics, would give new insights to its spreading history across Europe and allow the identification of relevant centres of dispersal.

Accounting for coypu in hunting laws would allow integrating it as a wildlife resource and harvesting coypu for its meat and fur. Meat of wild coypus was shown to be low in fat and cholesterol, while rich in proteins (Tulley et al. 2000) and the use of the coypu as a food source was common in Eastern Europe during the last century (Carter and Leonard 2002). Increasing the market value of the species would introduce an incentive for trappers and hunters and was shown to result in population decreases (Carter and Leonard 2002, Scheide 2013); however, it may also result in the wish to manage populations for permanent resource extraction.

Aside from direct control measures, another aspect of coypu management is the facilitation of winter survival and rise of reproduction rates by providing additional food sources. Wildlife feeding in or nearby settlements can induce rapid increases of coypu populations. Urban feeding sites and easily accessible agricultural areas may suffer from high coypu abundances (Walther et al. 2011, Scheide 2013). Managing urban populations is aggravated by the fact that the general public is often against lethal control methods of charismatic species, such as furry mammals (Walther et al. 2011, Jarić et al. 2020). Feeding bans, in combination with educational measures, such as awareness campaigns in such areas, are thus essential and have multiple positive benefits.

It is well-established knowledge that the coypu causes substantial economic and environmental damage when occurring at high densities. Although cool-temperate climates were believed to keep coypus at low densities, in many parts of Europe numbers have increased strongly during the last decades (Gosling and Baker 1989, Carter and Leonard 2002, Scheide 2013). Therefore, in Europe, the species was already subject to several national control campaigns (Carter and Leonard 2002), peaking in its inclusion as species of European Union concern (EU Commission 2014).

Our study shows that the coypu has, by far, not yet reached all potentially suitable regions in Europe and further highlights the importance of clarifying its response to increasing temperatures and arid conditions as they are likely to increasingly occur in the near future under climate change. However, one must consider the shortcomings of predictions that are made on the basis of opportunistic records from various sources and of differing data quality. Sampling effort differs spatiotemporally across the study area and, although we considered the violated assumption of an equilibrium for taxa undergoing an invasion process (Elith and Leathwick 2009, Václavík and Meentemeyer 2009, 2012) by depicting the uncertainty in predictions through a consensus approach, the outcome should be interpreted with caution (Pearce and Boyce 2006). Predictions of invasion processes should be regularly reassessed, ensuring that eventual changes in the species realised niche are captured (Václavík and Meentemeyer 2012). Nevertheless, SDM-based predictions of alien species’ distributions provide valuable tools to predict the next phase of invasion (Václavík and Meentemeyer 2012) and areas at high risk and can serve as the basis for further detailed analyses on regional or local scales, helping to better allocate resources for both surveys and management.