Global occurrence records for coypu were downloaded from the Global Biodiversity Information Facility (GBIF; gbif.org; March 4, 2011). The data were inspected, and records with a resolution greater than 30 minutes, our model resolution, were removed due to accuracy issues. We also removed presence locations in countries or states with a status of never established or extinct, retaining only those with a status of country of origin, escape or release, range expansion, or eradicated as defined by a global review of coypu distribution (Carter and Leonard 2002).

For the USA, we used monthly mean, minimum and maximum temperature data at a 4-km spatial resolution between 2003 and 2007 for our analyses (PRISM Group 2007). This time frame is biologically meaningful in that it matches the average lifespan of an individual, and is data-driven in that it matches the time frame of the sub-watershed scale (hydrologic unit code [HUC] 12) data used to validate the model in the Pacific Northwest, USA. Global environmental data were obtained from WorldClim (Hijmans et al. 2005). These data were averaged by month between 1950 and 2000 at a 30 arc second resolution. Thus, the national-scale climate data had a fine temporal resolution (monthly data) that matched some data collection whereas the global climate data had a finer spatial resolution with a coarse temporal resolution (50-year average).

We developed ecophysiological based models at the continental USA and the global scale, based on known physiological constraints on coypu. This species has known winter temperature tolerances that are thought to be the primary limiting factors on their distribution, at least in temperate regions (Gosling et al. 1983). Gosling et al. developed a population simulation model based on observed relationships between sequences of freezing days (defined as minimum temperature < 0 °C and maximum temperature < 5 °C) and body fat, litter frequency, litter size, and mortality. This model showed that a sequence of freezing days resulted in population declines due to adverse effects on the four measured characteristics, and Doncaster and Micol (1990) reported a 71% decrease in population density after canals were frozen for 20 consecutive days in France. In addition, coypu heavily depend on aquatic environments and are limited to environments within the transition zone between aquatic and upland environments (D’Adamo et al. 2000; Doncaster and Micol 1989; Guichón et al. 2003). We used this information on known coypu requirements (sensitivity to cold temperatures and need for aquatic environment) to define a model of habitat suitability rather than allowing a statistical model to detect relationships between habitat suitability and coypu presence.

For the continental USA we developed two different ecophysiological based models using monthly climate data from the PRISM data set at 4-km resolution; one using a five-year period (2003 to 2007) and another using a three-year period (2005 to 2007), hereafter referred to as US 5yr and US 3yr. We used two different time periods to assess the importance of inter-annual climatic variability on predicted distribution. We calculated the number of months within each time period that had a minimum temperature of less than 0 °C and a maximum temperature of less than 5 °C. Given the negative relationship between coypu populations and sequences of freezing days, we defined any month with average values meeting these criteria as unsuitable for coypu survival. To address the water limitation we developed a layer of arid locations by identifying locations in the USA with annual precipitation less than 250 mm, based on PRISM average annual precipitation from 2003 to 2007.

For the global ecophysiological model, we used WorldClim monthly data averaged from 1950 to 2000 at a 30 arc second resolution (~1 km), hereafter referred to as Global 50yr. Unsuitable environments were defined as locations with any month meeting the criteria of average minimum temperature less than 0 °C and average maximum temperature less than 5 °C. We again masked out arid regions, defined as areas with annual precipitation less than 250 mm based on the WorldClim average annual precipitation layer.

We used the VisTrails software (Freire et al. 2006) with the Software for Assisted Habitat Modeling (SAHM) package (Morisette et al. 2013) to develop correlative models of global coypu distribution using a Generalized Linear Model (GLM). We used GLMs because this technique creates simple models and has been recommended for model generalization (i.e., transferability to novel environment or time periods, Heikkinen et al. 2012). We generated background points using two different methods: a random generation of 10,000 locations within countries from which our data set had location records, referred to as GLM country, and a targeted background approach using location data for muskrats (Ondatra zibethicus), referred to as GLM targeted (Suppl. material 1: Figure 1). We downloaded muskrat data from GBIF and cleaned it by dropping records that had a spatial resolution greater than 30 min, removing fossil records and removing records from countries not known to have muskrats. The target background approach of using locations for similar species is recommended when using a presence-background method where the data are likely to have sampling bias (Phillips et al. 2009). Coypu and muskrats are sympatric species because both rodents are aquatic, are herbivores, are burrowers, and have similar global distributions (Ruys et al. 2011). By using a targeted background approach, biases in the presence locations are also assumed to occur in the background and thus cancel each other to some degree, similar to presence and absence data collected using the same methodology.

Figure 1.

USA state assessment of the five models. The assessment includes the number of USA states classified with at least some suitable (1) or no suitable (0) coypu (Myocastor coypus [Molina, 1782]) habitat for each coypu status class (never established/ extinct, present, no data, or eradicated) as defined by Carter and Leonard (2002) for each model of the five models.

The number of environmental variables from the global WorldClim data set used in the GLM was limited to six based on the known physiology of coypu and included mean diurnal range, maximum temperature of the warmest month, minimum temperature of the coldest month, annual precipitation, precipitation seasonality, and precipitation of warmest quarter. Environmental variables were reduced by removing one of each pair of highly correlated environmental variables (maximum of Spearman rank coefficient, Pearson’s product moment or Kendall tau rank; |r| > 0.7 following the recommendation of Dormann et al. (2013)) and biological knowledge of the species.

Using a threshold defined as maximizing sensitivity plus specificity divided by two, we created binary predictions of suitable and unsuitable habitat for the correlative models. Liu et al. (2013) recommended this threshold because it is transferable between methods that use presence-absence and presence-background. The binary predictions were then used to create equal-weight ensemble predictions of habitat suitability for coypu. We created two sets of ensemble models, three for the USA and one for the globe. The USA ensembles included an ensemble of the correlative models (GLM country and GLM targeted), another of the three ecophysiological based models (US 3yr, US 5yr, and global 50yr), and another of all five models. The global ensemble model was created using all three global models including the global 50yr, the GLM country, and GLM targeted models.

The models were evaluated using zonal statics at two scales; sub-watershed hydrologic unit code (HUC12) and the USA state boundaries. Standardized spatial surveys completed by on-the-ground fish and wildlife biologists for Oregon and Washington provided coypu density estimates at the HUC12 level and were used as an independent model validation (Sheffels 2013). Using the binary model predictions, we calculated zonal statistics using ArcGIS version 10.0 (ESRI 2011) for each HUC12. If any location within a HUC12 was classified as suitable by the model, the entire HUC12 was defined as suitable. HUC12 coypu density estimates were grouped within four density classes, >100, 11–100, 1–10, and 0 individuals, and the percent of HUC12 units that classified as suitable for each model were calculated for each density class. The models were also evaluated using zonal statistics identifying the number of USA states by coypu status (i.e., never established/extinct, present, no data and eradicated) identified by Carter and Leonard (2002). Again, if a state had any locations within it identified as suitable by the model, the state was defined as suitable, while a state was defined as unsuitable if it did not have any suitable habitat (i.e., no suitable locations within entire state).

We evaluated the global models using two additional methods. Similar to the state level evaluation, we used country level zonal statistics compared to the coypu status identified by Carter and Leonard (2002) for countries. Countries were classified into two coypu occurrence statuses: present (status of present or eradicated) and absent (status of never established or extinct). A country was classified as suitable habitat if any location within the country was predicted suitable based on the model. In addition, we used independent georeferenced locations as another evaluation metric. These independent records were compiled by searching the social media site ‘YouTube’ for the keywords: ‘bieberratte’ and ‘wasserratte’ (German), ‘beverrat’ (Dutch), ‘castorino’ (Italian), ‘coypu’ (British English, Spanish), ‘nutria’ (American English, Italian), ‘ragondin’ (French), Нутрия (Russian, Kyrgyzstani, Uzbekistani). The coypu in the videos had to be a naturally occurring population and the location of where the video was taken provided. Videos were examined to make certain that other species were not being misidentified as coypus. Videos where coypu were held as pets or in other confined situations such as fur farms, zoos or aquaria or for which location could not be determined were excluded. For videos not in English or French we used ‘Google Translate’ as an approximate translation tool to determine the circumstances and location. Since our focus was documenting the range of coypu, once presence was determined in a particular location we excluded videos from those regions on future searches. Finally, we used documents from national reports such as the ‘Red Book Data and Invasive Species Korea’ or personal communications from trusted researchers to verify regional presence. We used R version 2.15 (R Core Team, 2012) and the caret package (Kuhn 2013) to calculate sensitivity, specificity and percent correctly classified for each global model based on the model predictions and either the classification by Carter and Leonard (2002) or the independent locations.

We applied our ecophysiological based rule-set to future climate data. We obtained historic data from the Maurer data set (Maurer et al. 2007), which covers the USA at 1/8^th degree (~12km). These data were the reference data set used to downscale the global climate projections (GCMs) from the World Climate Research Programme’s Coupled Model Intercomparison Project phase 5(CMIP5) projections using the monthly bias correction spatial disaggregation (BCSD) technique (Reclamation 2013). These downscaled GCMs (listed in Suppl. material 1: Table 1) provided monthly projections of total precipitation and monthly average temperature. We quantified the amount of suitable habitat (i.e., the number of pixels meeting the criteria of monthly minimum temperature >0 °C or monthly maximum temperature >5 °C) for the Maurer dataset (1950 to 2013) and the 12 downscaled GCMs available (1950 to 2013 for historic comparison; 2014 to 2100 as climate forecasts). We calculated the forecasts for all four representative concentration pathways (RCPs), which describe four different greenhouse gas concentration trajectories. We also generated maps of suitable habitat based on the Maurer dataset for 2003 to 2013 and an ensemble of downscaled GCMs for 2006 to 2016 to assess how well the GCMs performed currently when compared to the observations on which they were calibrated. Because the forecasted GCMs began in 2006, we were unable to have complete overlap in the decades for comparison. We applied our criteria for suitable habitat to the GCMs for the period 2040 to 2050 to assess how coypu distribution may change in the future. For this forecast we only used the 4.5 RCP, as RCPs do not begin to diverge significantly until after mid-century.

For the ecophysiological based models, we produced layers with the number of months for each cell that did not meet the required temperature criteria. For the US 5yr model, the number of months with unsuitable temperature conditions ranged from 0 to 41, while for the US 3yr model the maximum number of unsuitable months was 28.

The GLM country model retained all six environmental variables in model fitting, while the GLM targeted model dropped average annual precipitation and mean diurnal range. Average minimum temperature of the coldest month was the most important predictor in both models, with a logistic shape where suitability began to steeply increase from zero around -10 °C and climbing to 1.3 °C before reaching an index value defined as suitable. Both models retained temperature of the warmest month, with a generally positive relationship when considered with other variables. However, a function considering that predictor alone revealed a hump shaped relationship. Internal cross validation produced good assessment metrics for both models. The GLM country model had a cross-validation area under the receiver operating characteristic curve (AUC) of 0.94 and a true skill statistic (TSS) of 0.76, while the GLM targeted model had a cross-validation AUC value of 0.91 and TSS of 0.70. To produce binary maps, the GLM country threshold was 0.14 and the GLM targeted threshold was 0.44.

All models performed well when compared to the HUC12 coypu density data (Table 1). All HUC12 sub-watersheds with a coypu density greater than 100 individuals per sub-watershed were predicted by all models to contain suitable habitat. The percentage of coypu density HUC12 classes of 11-100 and 1-10 predicted as suitable by the models were also high (> 87%), while HUC12 areas with a reported density of 0 had a much lower percentage predicted as suitable (Table 1). Model predictions compared to USA state-level classifications showed the models were better at correctly identifying states with established coypu populations than in predicting states that had populations that are now extinct or states where coypu have never been established (Fig. 1). Only the US 3yr and US 5yr models predicted a state classified as having an established coypu population as unsuitable (state of Delaware), while the other three models correctly classified all states with coypu status as present. For states with no data on population status, three to 12 of them were predicted to contain suitable habitat.

The ecophysiological based ensemble model results show greatest agreement in suitability in the southeastern USA from Texas to North Carolina and along the Pacific coast from Washington to southern California (Fig. 2a). Currently, coypu are not in California. We hypothesize this is due to a geographic barrier to their expansion south from Oregon. The waterways containing coypu in Southwestern Oregon are not hydrologically connected to the ones in Northern California and mountain ranges separate the two. The greatest area of model discrepancy was along the border between Tennessee and Arkansas, where the US 3yr model predicted a more northern distribution limit compared to the US 5yr and global 50yr models (Fig. 2a). These models agreed on suitable/ unsuitable classification 92% of the time. The GLM targeted model predicted a much greater amount of suitable habitat across the southern USA than the GLM country model (Fig. 2b). The ensemble of all five models for the USA revealed that the addition of the correlative models to the ecophysiological based models resulted in a more restricted distribution of agreement (Fig. 2c). Model agreement on habitat suitability was more confined to east Texas, Louisiana, Mississippi, and Alabama, mainly following the restricted distribution of the GLM country model. On the Pacific coast, model agreement was restricted to areas in Washington, Oregon and the northern half of California.

Figure 2.

Model predictions for Myocastor coypus [Molina, 1782] for a an ensemble of US 3yr, US 5yr, and global 50yr b an ensemble of GLM country and GLM targeted c an ensemble of all five models d number of months classified as unsuitable using the Maurer observed climate data for 2001 to 2010 e the number of GCMs defining each pixel as suitable (ensemble of the 29 binary downscaled GCMs using the Maurer dataset as the reference) f ensemble of the 31 downscaled GCMs average from 2040 to 2050. All maps are overlaid with USA state population status according to Carter and Leonard (2002). Unsuitable habitat is defined as areas where no models predicted the area as suitable, while suitable areas are defined according to which model(s) predicted suitable habitat. Maps are in Albers Equal Area projection.

At the global scale, global 50yr, GLM targeted, and GLM country models had varying levels of performance when compared to country level classification by Carter and Leonard (2002; Table 2). For country level comparison, all models had high sensitivity values (>0.85) and low specificity values (<0.25) with the number of countries correctly classified > 70% (Table 2). When comparing global model performance to independent coypu locations, sensitivity values decreased to 0.73 for all models and specificity values and percent correctly classified increased (Table 2).

An ensemble of the three global models had high agreement in coypu suitability for regions with established invasions (Fig. 3), although all three models only agreed on classification as suitable or unsuitable across 59% of the globe. The GLM targeted and global 50yr models were more similar with 81% agreement between predictions. The GLM targeted model predicted much more of the earth’s surface as suitable (59%) compared to the other models (GLM country = 32%; global 50yr = 45%). All models predicted suitable coypu habitat in Western Europe and portions of the USA where coypu populations are known to exist – and where occurrence data were available to fit the models. At least a portion of all countries reporting coypu as native were predicted as suitable by all models. Evaluation of the models in the native range is difficult, however, as range maps that do exist do not provide information on how they were derived, and the rigor with which they were created is questionable.

Figure 3.

Global predictions of habitat suitability for coypu (Myocastor coypus [Molina, 1782]). Ensemble predictions using three models at the global scale including an ecophysiological based model based on average monthly climate data from WorldClim, a correlative model using country background and a correlative model using a taxonomically targeted background approach. Maps are in Mollweide projection.

The GLM country model predicted the least amount of tropical areas as suitable, with the GLM targeted model and the global 50yr model being more similar. However, many of these areas had novel environmental conditions. In dry areas such as North Africa, however, the global 50yr model did not predict suitable habitat due to the added arid region mask. The GLM country model excluded some of these dry areas, while the GLM targeted model included almost all of them.

Model evaluations from HUC coypu density for the northwestern USA show very little difference between ecophysiological based and correlative models (Table 1). The same was true for evaluations at the state level. The greatest difference was for the ‘no data’ category where the GLM country model predicted suitable habitat in only ten states while the other models predicted suitable habitat in three to six states (Fig. 1). For the global models, predictions between the correlative models and the ecophysiological based model were again very similar for both country level evaluations and independent location evaluations.

There is substantial variation in potential future climate from year to year as well as between GCMs and RCPs (Fig. 4). RCPs do not vary greatly until after mid-century. Comparing predictions of suitability between the modeled climate for current conditions (2006-2016) and the observed climate on which it is based (Maurer 2001-2010) revealed some discrepancies, with a mean disagreement in predictions (over or under prediction) of 6.1% of the area of USA (range of 3.2 to 10.6%). The Maurer dataset defined 34% of the USA as suitable (Fig. 2d). Extremes among the GCMs ranged from giss-e2-r-cc predicting 27% to canesm2 predicting 40% of the USA as suitable during the 2006-2016 period, though both had higher than average levels of disagreement (7.4% and 10.6%, respectively). The average amount of suitable area (32.4%) was comparable to the Maurer reference dataset. The range was similar for 2040 to 2050, with a minimum of 27.5% (ACCESS1-0), a maximum of 44.1% (CSIRO-MK3-6.0) and an average of 35.5%. Only three of the 29 GCMs predicted a decrease in suitable habitat (-3.7%, -0.2%, and -0.2%) and the maximum predicted increase in suitable habitat was 9.4% by HADGEM2-AO. The average increase in suitable habitat was an additional 3.1% of the USA. There was more discrepancy among GCM predictions for the eastern USA than the western USA (Fig. 2e and f).

Figure 4.

Amount of suitable habitat for coypu by year starting in 1950 and extending to 2100. Amount of suitable habitat is defined as thousands of km² within the continental USA without any months where average minimum temperature was <0 °C while average maximum temperature was also <5 °C. The solid black line from 1950 to 2013 is the Maurer observed dataset, the historical data is the 12 General Circulation Models (GCMs) calibrated between 1950 and 2013 using the Maurer dataset, and the projected climate by the GCMs with the average amount of predicted suitable habitat (solid line) and variation in predicted suitable habitat (solid colored area) for the four different representative concentration pathways (RCPs) describing possible climate futures by the GCMs. The solid vertical bars indicate the time periods for which we created geographic maps of predicted suitable habitat.