University of Amsterdam
BACHELOR THESIS
Nest-site preference of the muskrat Ondatra zibethicus in North Holland
A thesis submitted in the fulfillment of the requirements of
the bachelor’s degree of Future Planet Studies
June 29, 2018
Author:T. E. van Noppen
Supervisors: E. E. van Loon & C. Ootes
Abstract
The muskrat Ondatra zibethicus is considered a pest species in the Netherlands due to its digging activities which result in subsidence of dikes and unstable flood defenses. Pest control programs have been implemented by the Dutch Water Authorities in order to minimize the muskrat population. This research used an ecological-niche factor analysis to examine which locations are primarily preferred by the muskrat for nest construction in the district ‘Schermer’. Therefore, maps containing information about the elevation and soil and vegetation type were transformed into eight different eco-geographical variables which have been compared to presence data of nest-sites. The results indicate that the locations used for nest construction differ from the entire research area. In addition, the analysis shows that the muskrat prefers a sloping river embankment, pointing south-east, that consists of clay or fine sand and is covered by reed vegetation. However, only 21 reliable presence points for nest-sites about the research area could be used, thus further research is essential to make the results more robust.
Table of Contents
Introduction 4 Theoretical framework 5 Method s 7 Results 11 Discussion 15 Conclusion 17 Cited literature 18 Appendix 1. GIS-commands 20Appendix 2. Used maps 24
Appendix 3. Created maps 27
Appendix 4. Script RStudio 29
Introduction
The muskrat Ondatra zibethicus is a semi-aquatic rodent and is considered an invasive species in Europe. The animal originates from North-America and was introduced into central Europe in 1905 as a furbearer. Subsequently, through natural dispersal and further introductions in Europe the animal rapidly colonized large areas of the continent (Ecke et al., 2013; Skyriene & Paulauskas, 2012). In the Netherlands the first muskrat was discovered in 1941 (Heidinga, 2006). Currently, the animal is found throughout the country, except from the Wadden Islands Texel and Vlieland (Heidinga, 2006).
The muskrat is one of the animals with the highest potential economic impact of all alien mammals in Europe, which is mainly due to its digging activities (Ecke et al., 2013). For the construction of nests, the muskrat digs holes into earthen embankments, causing subsidence of dikes and unstable flood defenses, which results in dangerous situations for humans and cattle. Furthermore, the animal damages cultivated farmland in rural districts due to its gnaw activity (Heidinga, 2006). Hence, a year-round pest control program has been implemented by the Dutch Water Authorities to minimize the muskrat population and to keep damage and required repair costs below a publicly acceptable level (Bos & Ydenberg, 2011). This program led to a significant decline of the number of trapped muskrats from more than 400.000 in 2005 to 100.000 in 2016 (Unie van Waterschappen, 2016). Nevertheless, due to the high reproduction rate and the inflow of muskrats originating from surrounding countries, it is very difficult to eradicate the entire population in the Netherlands (Heidinga, 2006).
Controlling the muskrat population has been done by killing or trapping the animals during the whole year over the entire country, with one exception (the National park: the Oostvaardersplassen). This has been quite intensive, for instance in 2007 on average 14 man-hours per km2 were made and about 35 million euro was spent (Bos & Ydenberg, 2011).
This study is an attempt to examine the possibility of predicting nest location preference of the muskrat in the district ‘Schermer’ in the province North Holland. Understanding which areas are preferred by this animal is crucial in order to make pest control management more efficient, hence less expensive while decreasing animal suffering. The following main question has been used for this research:
To what extent is it possible to predict nest-site preference of the muskrat for the district ‘Schermer’ in North Holland?
In literature different types of analysis have been used to determine habitat preference of various species. Among these is the ecological-niche factor analysis (ENFA), which is considered particularly advantageous for this type of analysis, since it requires no absence data but presence data about the species of interest, which is compared to environmental conditions (Galparsoro et al., 2009). According to Hirzel et al. (2002), absence data of a species is very hard to obtain accurately. Hence, the ENFA has been selected as the optimal method for this research. Hitherto, no predictions with the use of ENFA about habitat suitability of the muskrat in the Netherlands have been made in scientific literature. However, this analysis could offer a significant aid in better predicting the preferred habitat of the muskrat.
To perform the ENFA, firstly, literature research has been carried out to gain insight in the environmental features that could be relevant for the analysis. In addition, the muskrat management organization gave insight into their control policy and demonstrated typical sites containing burrows. The information gathered from this research is presented in the next chapter.
Theoretical framework
For the theoretical framework, some background information about the muskrat will be given, so that it is clear what characteristics this animal has and how this relates to the environment the animal prefers. Subsequently, certain typical habitat features of muskrats according to literature will be explained. Finally, a brief description about the current approach of the pest-control program in the Netherlands will be given.
The muskrat
Muskrats are medium-size rodents with a semi-aquatic lifestyle and a strong social organization (Bos & Ydenberg, 2011). The exact size of the muskrat population in the Netherlands remains unclear. Nonetheless, it is certain that the muskrat population shows strong seasonal fluctuations, which is the result of the high reproduction rate of the animal, with on average three litters of about six young per year, and the high mortality rate, especially during fall and winter (Doude van Troostwijk, 1978; Van Loon et al. 2016).
In general, it is stated that the muskrats do not have any natural predators in the Netherlands, hence they have found a so-called ecological niche (Heidinga, 2006). Nevertheless, certain muskrats are predated by the fox (Vulpes vulpes), the American mink (Neovison vison), the Eurasian otter (Lutra lutra) and the European polecat (Mustela putorius). However, Heidinga (2006) states that the impact of these predators on the muskrat population in the Netherlands is negligible.
The muskrat is a rather strict herbivore, its diet includes more than 50 natural plant species, and several cultivated species, such as grain, maize, sugar beet, carrot and endive (Heidinga, 2006). During the summer the animal prefers to eat Typha (Typha latifolia), while in the winter it consumes mainly the roots of reed (Phragmites communis and Phragmites australis) and in certain situations when plant food is scarce, it feeds on freshwater mussels (Anodonta cygnea) (Ecke et al., 2013).
Plant materials are not only used as food, but are also utilized for house building (Danell, 1979). Depending on the environment, the muskrat constructs either a lodge or digs burrows in riverbanks (Heidinga, 2006). Lodges are conical or irregular mounds above the water level of about 1 – 2 meter in diameter with internal chambers (Connors et al., 2000). Burrows are tunnels of 13 – 15 cm in diameter and up to 13 meters long (Connors et al., 2000).
Preferred habitat features described in literature
Several characteristics about the preferred habitat of the muskrat have been discussed in literature. The first and most significant is the presence of water, as the muskrat is a semi-aquatic animal. Hence, Doude van Troostwijk (1978) states that the Netherlands offers a well-suited environment for the muskrat with its vast network of waterways. Secondly, Skyriene & Paulauskas (2012) indicate that the size of the muskrat population largely depends on the composition and the amount of available food. Furthermore, it is argued that muskrats are sensitive to changes in water level and population density (Doude van Troostwijk, 1978). In addition, water quality could also be an influential factor. Skyriene & Paulauskas (2012) state that habitat preference of the muskrat depends on the water pH, salinity and the amount of dissolved oxygen. However, muskrats have also been found in strongly polluted and brackish water (Doude van Troostwijk, 1978).
Allen & Hoffman (1984) state that muskrats prefer still or low velocity water. Which is in accordance with Heidinga (2006), who argues that muskrats prefer clear, stagnant water with a depth of
0.46 to 1.20 meter. Besides, it is argued that the most suitable river bank for burrow construction is a bank with a slope of 10 degrees or more and a minimum height of 0.5 meter (Allen & Hoffman, 1984).
Ecke, Henry and Danell (2013) were able to develop a model that could successfully predict muskrat occurrence in Swedish lakes with the use of topography and vegetation data. According to this model, the areas that were mainly suitable for muskrats included lakes with extensive areas of meadows rich in herbs and lakeshore meadows rich in halophyte vegetation.
Pest control management in the Netherlands
The main objective of pest management in the Netherlands has been to minimize the population (Bos & Ydenberg, 2011). During the whole year muskrats are captured and killed, which is done by the use of about ten different types of traps, and incidental by the use of a rifle. Other methods such as chemical or biological methods are not applied in the Netherlands (Van Vliet & Lengkeek, 2007). Trapping can be active or passive (Van Vliet & Lengkeek, 2007). Active trapping is a method whereby research is done about the exact location of a burrow, subsequently traps are placed at the entrances of the burrow by the use of a conibear trap (figure 1). Passive trapping indicates the installation of traps in watercourses to capture migrating muskrats.
Since 1987 all trapped muskrats in the Netherlands have been registered in an online data set, de vangstnet registratie (Heidinga,
2006). This data is used to examine the development of the number of captured animals per hour, which is used to obtain information about the efficiency of the muskrat control management (Heidinga, 2006). The total number of trapped muskrats is visualized in figure 2. The data appears to suggest that the pest-control program has resulted in a significant decline in the total muskrat population in the Netherlands.
x 1.000
Methods
The input of the ecological-niche factor analysis consists of presence data of muskrat nests and various eco-geographical variables describing the surroundings of the research area. Therefore, maps about environmental characteristics of the study area have been collected, which then were analyzed using data on muskrat catchments in the study area. Firstly, the ENFA will be further elucidated.
Ecological-niche factor analysis
The majority of studies about habitat suitability of a certain species have been conducted using logistic regression, which contains a dependent variable in the form of presence/absence data (Xuezhi, 2008). However, research has shown that it is very difficult to accurately determine the absence data, which might result in the underestimation of habitat preferences and a higher probability of type II errors (Ecke et al., 2013). Thus, for this research ENFA will be used, whereby only presence data is required. ENFA compares, by using eco-geographical variables (EGVs), the locations where the focal species has been observed to a reference set describing the whole study area (Xuezhi, 2008). EGVs may represent topographical features (e.g., altitude, slope), ecological data (e.g., nitrate concentration, frequency of forests), or human structures (e.g. distance to the nearest town, road density) (Hirzel et al., 2002). ENFA summarizes all predictors from the EGVs into a few uncorrelated factors, and then builds a habitat suitability map based on eigenvectors and eigenvalues (Xuezhi, 2008). ENFA builds on three conceptions: Marginality,
Specialization,
Tolerance,
where mG represents the mean of the EGV in the whole study area, σG is the standard deviation of the global distribution, mS is the mean of the EGV of the species distributions and σS is the standard deviation of the focal species distribution. The coefficient weighting σG (1.96) ensures that marginality will be most often between zero and one (Hirzel et al., 2002). Namely, if the global distribution is normal, the marginality of a random chosen cell has only 5% chance of exceeding unity (Hirzel et al., 2002). Marginality describes the ecological distance between the species optimum and the mean habitat in the whole area (Xuezhi, 2008). A positive value of marginality indicates that the species prefers values higher than the mean in the entire studied area, while a negative value indicates the preference of lower-than-mean values. Specialization has been defined as the ratio of the ecological variance in the mean habitat compared to that of the investigated species (Xuezhi, 2008). Specialization varies from one to infinite, but the tolerance value, which indicates the inverse of the specialization, takes a value between zero and one.
A value close to zero indicates that the species are distributed in a narrow space, while a value close to 1 indicates that the distribution of the species is large (Xuezhi, 2008).
Study area
This research has been performed in North Holland. where the number of trapped muskrats per hour is relative low (LCCM, 2007), suggesting a low population density in this area. Measuring muskrat density in the research area is very difficult and will not be included in this research, however population density could be of influence on the nest-site preference of the muskrat. Therefore, a low-density area has been opted as the optimum research area of this analysis. In addition, to narrow down the
scope of the project a sub-region of North Holland, the district ‘Schermer’, has been selected for this research. Which is a district consisting of 64.87 km2 and is located in the center of North Holland (figure
3).
After the selection of the study area, a mask layer has been prepared in ArcMap (version10.4.1), which could be used to extract only the relevant data from the various maps. The mask layer consisted of cells of containing water, or cells located in the surrounding of 25m of a water body, since nest-sites are usually encountered in or nearby water, the resulting map is visualized in figure 4. For the step-by-step process in ArcMap, including all the spatial analysis tools that have been used for this research, is referred to appendix 1.
The cell-size of the used for the mask layer was 5x5m, which is the cell-size that has been for all the maps that were used for the analysis. For analyses about habitat suitability, the spatial resolution determines largely the accuracy and reliability of the obtained results (Chefaoui et al., 2005). The used cell size of 5x5m made the size of the datasets, in terms of size, appropriate for the ENFA, while keeping the variation of the variables for each grid-cell as small as possible. After extracting the relevant data, all the maps have been imported into RStudio (version 1.1.4), where the ENFA has been carried out.
Figure 4: prepared mask layer. Blue area indicates water and areas in the surrounding of 25m of a water body.
Data on muskrat occurrence
For this research, the “presence” data of muskrats has been collected using the dataset the ‘vangstnet
registratie’ which has been obtained from the Dutch Water Authorities of North Holland. This dataset
contained data about all the trapped muskrats between 23-05-2014 until 14-06-2018 and provided information about the exact time and the location of the captured muskrats, which is visualized in figure
5.
The dataset has been imported into ArcMap in order to construct a binary map of the nest locations of the muskrat. Therefore, a selection has been made for the muskrats that were trapped with a conibear trap, which is the type of trap that is used to catch the animal at the entrance or in the surroundings of a burrow. Furthermore, a selection has been made for the catchments that where performed inside the research area. Consequently, the obtained data points were converted into a binary map, with ‘1’ indicating the areas where a nest has been found, and ‘0’ indicating the cells where no nests were encountered. The resulting layer was imported into RStudio, so that it could be used as presence data for the ENFA.
Eco-geographical variables
Literature research and interviews with staff members and the coordinator of the muskrat control organization in North Holland revealed a variety of eco-geographical variables that could be relevant for this research. However, the final selection of the eco-geographical variables that has been used for this research was mainly based on the available data about the investigated area. The used data has been obtained from an online database of the University of Amsterdam and from the Dutch Water Authorities
of North Holland (Hoogheemraadschap Hollands Noorderkwartier) The EGV’s that eventually have been used for the ENFA are displayed in table 1. For an overview of the used maps is referred to appendix 2 and 3.
Subsequently, the three maps have been imported into the program ArcMap to transform them into raster-layers, hence making them appropriate for analyzation in RStudio. Therefore, the digital elevation map has been converted into a slope map, representing the rate of change of elevation for each grid cell, which was used to determine the slope of the river embankments in the research area. Additionally, an aspect map has been prepared with the use of the digital elevation map. The aspect map has been used to establish the direction of the slope of the river embankments, which could be used to determine whether nest site preference would be related to the amount of sun a river embankment receives.
The land cover and soil map consisted of a variety of polygon features. Both maps were converted into raster layers with a cell-size of 5x5m and only relevant data was subtracted from the layers with the use of the prepared mask layer. This mask layer has also been used to extract data from the aspect- and slope layer, hence all four maps contained the same cell-size and number of cells. Afterwards, the maps were imported into RStudio. With the aid of RStudio, six derived EGV’s were obtained from the land cover and soil map, namely four different types of soil maps, that were used to examine the influence of soil type on the nest site preference. And two different types of vegetation maps, which have been used to investigate the influence of reed and tree vegetation. For the R-script that has been used for this research is referred to appendix 4. In the next chapter the resulting eight map are visualized and the outcome of the ENFA will be discussed.
Type of data Derived EGV Description
Digital Elevation Model (source: Algemeen
hoogtebestand Nederland (ANH2), available at the geoportal of the University of Amsterdam)
Slope map Rate of change of elevation for grid cell of a digital elevation model (DEM). First derivative of a DEM. Values range between 0-90 degrees.
Aspect map Indicates downslope direction of the maximum rate of change in value from each cell to its neighbors. Measures clockwise in degrees from 0 (pointing northwards) to 360 degrees (again direction towards the north).
Land cover
(source: Hoogheemraadschap
Hollands Noorderkwartier)
Reed vegetation Presence or absence of reed vegetation. Indicated with a value of ‘1’ for presence and ‘0’ for absence. Tree vegetation Presence or absence of tree vegetation. Indicated
with a value of ‘1’ for presence and ‘0’ for absence. Soil map
(source Hoogheemraadschap
Hollands Noorderkwartier)
Sandy soils Presence or absence of sandy soils. Indicated with a value of ‘1’ for presence and ‘0’ for absence. Bog soils Presence or absence of peaty soils. Indicated with
Fine sand soils (‘zavel’)
Presence or absence of soils containing between 8-25% of small particles (2 µm). Indicated with a value of ‘1’ for presence and ‘0’ for absence. Clay soils Presence or absence of clay soil, which is defined
as a soil containing >25% small particles (2 µm). Indicated with a value of ‘1’ for presence and ‘0’ for absence.
Table 1: Eco-geographical variables used for analysis
Results
In the period of 23-05-2014 until 14-06-2018 a total of 13.821 muskrats have been trapped in the province North Holland. The majority of these animals, in total 9.314, were captured with the use of a conibear trap. After selecting on the district ‘Schermer’ and omitting the grid cells that contained no data, only 21 reliable presence points could be obtained. Subsequently, the data was imported into RStudio to visualize the obtained maps and to perform the ENFA.
Visualization eco-geographical variables
For the EGV’s several kernel density plots have been constructed in order to obtain an overview of the distribution of the eight variables. Additionally, a selection has been made for those grid cells that contained one of the 21 samples, which have also been visualized with the use of a kernel density plot. This resulted in eight different graphs indicating both the density estimation of the entire research area and the density estimation of the areas that contained burrows. In figure 6 the plots are displayed, indicating also the various bandwidths that were applied. For a better visualization of the plots is referred to appendix 5.
Regarding the eight kernel density plots, several initial conclusions can be drawn. Firstly, the plots reveal that the reed and tree density is very low in the entire research area and for the grids that contained nest-sites. A similar distribution is visible in the density plot of sandy soils, indicating a very low proportion of grid cells consisting of sandy soils. The fine sand and clay plots are more evenly distributed for the both the reference area and for the 21 selected grid cells. The bog soil shows an even distribution for the reference area, however the distribution changes if a kernel density is plotted for the 21 grid cells. Furthermore, the slope plot shows that the investigated area is relatively flat, since the plot has a very right-skewed distribution. However, certain peaks above 10 degrees are visible in the nest-site plot, hence it could be argued that the muskrat prefers constructing a burrow in a steeper riverbank. Therefore, further analysis is necessary to examine whether the slope variable is linked to the nest-site preference of the muskrat.
Results ecological-niche factor analysis
The ENFA in RStudio resulted in an overall marginality of 1.33. According to Chefaoui et al. (2005), when the overall marginality obtains a value higher than 0.65, it can be stated that there is a difference between sample sites and the reference area.
In addition, the analysis returned probability plots for each of the EGV’s, which are similar to the kernel density plots, however, the probability plots of the ENFA indicate also the mean values of the EGV for the entire research and of the values of the 21 sample sites (figure 7).
Furthermore, the output of ENFA contained the values for marginality and suitability for each of the EGV. These are represented in table 2. The marginality scores indicate how much the values of each EGV of the nest locations differed from those of the entire research area (0 indicating not much difference from the mean and -1 and 1 indicating the preference of respectively lower and higher values than the mean of the entire region). The marginality coefficients show that the nest-sites of muskrats are essentially linked to aspect, slope, reed and the soil types: bog and fine sand (aspect = -0.61, slope = 0.28, reed = 0.24, bog = -0.46, fine sand = 0.75). Consequently, it could be argued that the muskrat prefers sloping river embankments covered with reed vegetation and with soils containing fine sand for the construction of a burrow. On the other hand, the animal seems to prefer areas containing no bog soils. Additionally, considering the mean values for aspect in figure 7, it could be stated that the muskrat prefers river embankments pointing south-eastward.
The values for specialization of the nest-site preference is mainly conditioned by the aspect (0.52) and clay (0.48) and bog (-0.65) soils, indicating that the variation of the EGV-data of these three variables at the nest-sites differed from the variation of the EGV-values of the reference area, consequently indicating that the nest locations show some specialization compared to the entire research area.
EGV Marginality Specialization
Aspect -0.60876792 0.51714993
Slope 0.27787947 -0.18080070
Clay (Soil) -0.19063140 0.48104360
Bog (Soil.1) -0.46425526 -0.64612057
Sand (Soil.2) -0.09928665 -0.12314857
Fine sand (Soil.3) 0.75024967 0.17068869
Tree (Begr_ter) -0.04589113 -0.05692030
Reed (Begr_ter.1) 0.23664132 0.05875541
Table 2: Marginality and Specialization of EGV's
To visualize the outcome of the ENFA a biplot has been produced, which is represented in figure 8. Here, the x-axis describes the amount of marginality and the y-axis the amount of specialization. The values displayed in table 2 have been used as the coordinates for the arrows, hence indicating the marginality and specialization of the 21 sample sites. The light and dark areas correspond to the minimum convex polygon enclosing all the projections of the available and used points respectively (Basille et al., 2008). The white dot on the y-axis corresponds to the centroid of the used nest locations of the muskrat, which is located slightly to the right of the axis, showing that the preferred nest location of the muskrat is somewhat different from the mean available conditions in the entire research area.
Discussion
The environmental niche factor analysis revealed that the muskrat tends to construct a burrow in slightly sloping river banks containing reed vegetation and soils consisting of fine sand. However certain remarks about the procedure of this research and the ENFA must be made regarding the results.
Obtained data
Firstly, several remarks could be made regarding the obtained catchment data. According to Hirzel et al. (2002), the presence data should consist of unbiased samples, representing the actual distribution of the focal species. For this research, the presence data was based on a dataset provided by the muskrat organization of North Holland. Uncertainty exist about the accuracy of this dataset. The obtained muskrat burrow distribution could have been slightly biased. For instance, due to the fact that certain areas might have been examined more frequently than others, since certain areas are easier to attain and perform field surveys than others (Galparsoro et al., 2009).
Secondly, it could be argued that the obtained data about the trapped muskrats is sensitive to small errors. The data is gathered by various muskrat fighters, hence mistakes about the exact location could easily enter the dataset, for instance due to the fact that the registration of some trapped animals occurred one or two days after the muskrat had been captured. Furthermore, the final sample size that was used for this research contained 21 samples, which is not very high. Further analysis could be performed in different districts with larger sample sizes. Additionally, other eco-geographical variables could be included to obtain more robust results, for instance about water depth or the size of the water body. Especially variables that show significant variation over the reference area would be useful. For example, in this research the amount of tree vegetation in the ‘Schermer’ was very low. It could be argued that another EGV would have been more relevant for this research. However, due to the small amount of available data about the research area these eight EGV’s where selected as the most valuable.
In addition, for the EGV’s no research has been performed about a possible correlation between certain EGV’s. All the variables have been used for the analysis, however, certain variables might have shown correlation and consequently could have been omitted of the ENFA, since correlating variables explain a similar percentage of variance (Galparsoro et al., 2009).
Preparation maps
For the preparation of the data, both the size of the study area and the grid-size were of great importance for the accuracy of the results. It could be stated that a larger reference area would probably have resulted in higher values for the marginality and suitability. The ENFA characterizes ecological niches relative to an in advance specified reference area, hence marginality and specialization are depended on the geographic limits of the research area (Hirzel et al., 2002).
Another point to emphasize is the grid size that has been used, in this case a cell size of 5x5m was used. However, the Digital Elevation Model was based on a cell-size of 0.5x0.5m. Several attempts using various cell-sizes have been performed in RStudio to determine that a minimum size of 5x5m was required so that the entire analysis could be performed. It could be argued that a smaller cell-size would obtain more accurate results. Nevertheless, the polygon features of the land use map and the soil map demonstrated a lower accuracy than the Digital Elevation Model. Therefore, the 5x5m turned out to be the most appropriate for this research.
A final remark about the preparation of the data must be made concerning the extensions of the various raster layers. After the datasets were imported into RStudio the layers contained slightly different extents. Several attempts with various tools in ArcMap did not obtain the desired similar extents, hence, some modifications in RStudio had to be made to the land cover and soil map, resulting in four maps that
could be used for the analysis. However, this could not avoid the fact that each grid cell of the land cover and soil map differed one meter from the grids of the aspect and slope map.
Ecological-niche factor analysis
The ecological-niche factor analysis has been very useful for this research, sine the analysis only requires presence data which decreases the chance of using ‘false’ absence data, due to the animal’s ability to disperse or hide during field surveys (Galparsoro et al., 2008). In addition, the ENFA aided in comparing various environmental features to the presence data of nest-sites. In this respect, the ENFA differs fundamentally from other analyses for instance first-order regressions and discriminant functions where relationships are presumed to be linear and monotonic (Hirzel et al., 2002).
However, certain limitations regarding the ENFA are discussed in scientific literature and are worthwhile to mention. Firstly, the results could be obviously affected by another feature that has not been included in this research (Hirzel et al., 2002). For instance, it might appear that reed vegetation is a factor that influences the preference of nest location of the muskrat, however, in reality it could turn out that it is not the reed the animal is looking for, but another type vegetation which is often located in the surroundings of reed plants. Nevertheless, it could be stated that the results provide useful insights about preferential conditions (Hirzel et al., 2002).
In general, it could be stated that habitat models are important tools that can be used to better comprehend habitat suitability of a certain species. Nonetheless, one must not forget that these models are no exact representation of reality, thus the interpretation of these kind of models must be done carefully (Galparsoro et al., 2008). However, they can provide valuable information for resource management and habitat conservation. Several results of this analysis are consistent with the variables that were mentioned by the muskrat fighters, for instance about the presence of reed which is used by the muskrat as feed. The results of this research could be used by the muskrat fighters of North Holland and other muskrat control organizations in the Netherlands. However, further investigation about nest preference or habitat suitability of the muskrat is essential, since the interviews with the muskrat fighters revealed a variety of features that could be relevant to examine, for instance about the type of vegetation that is usually encountered in the neighborhood of a burrow. Hitherto, not much scientific research has been done about these features. In addition, the results obtained from this study could be verified by performing a similar type of analysis in another district of the Netherlands.
Conclusion
For this analysis eight different eco-geographical variables have been compared to the presence data of nest locations of the muskrat. A variety of environmental features about habitat preference of the muskrat were encountered in scientific literature and have been taken into account while selecting the EGV’s that would be used for this research. The variables that were used are: slope, aspect, reed and tree vegetation and sand, bog, clay and fine sand soils. Subsequently, the appropriate cell-size and study area have been selected and the presence data has been prepared by extracting the muskrat catchments that were performed with a conibear trap in the district ‘Schermer’.
To obtain insight about the distribution of the data, Kernel density plots for all the variables where performed, revealing that the elevation of the study area did not show much variation in elevation, furthermore, the plots indicated that the area contained mainly clay and fine sand soils, and, according to the land use map, tree and reed vegetation were low in the research area.
After transforming the data into raster layers, while using a mask layer including only the areas containing water or those located in the surroundings of water, the ENFA has been applied to the eight maps and the presence data, which resulted in an overall marginality of 1.33 indicating that the 21 sample locations differed from the reference area. Furthermore, the ENFA revealed that the nest location of the muskrat is primarily linked to aspect, slope, reed vegetation and the soil types bog and fine sand.
To conclude, it is possible to predict nest-site location for the district ‘Schermer’, hence making pest control management in this district more efficient. However, due to the small sample size and the low number of EGV’s that were used for this analysis, the results are not very robust. Therefore, further research is essential to make the results more robust.
Cited literature
Allen, A. W. & Hoffman, R. D. (1984). Habitat suitability index models: muskrat. National Wetlands
Research Center.
Basille, M., Calenge, C., Marboutin, E., Andersen, R. & Gaillard, J.M. (2008). Assessing habitat selection using multivariate statistics: some refinements of the ecological-niche factor analysis. Ecological
modelling, 211, 233-240.
Bos, D. & Ydenberg, R. (2011). Evaluation of alternative management strategies of muskrat Ondatra
zibethicus population control using a population model. Wildlife Biology, 17(2), 143-155.
Connors, L. M., Kiviat, E., Groffman, P.M. & Ostfeld, R. S. (2000). Muskrat (Ondatra zibethicus) disturbance to vegetation and potential net nitrogen mineralization and nitrification rates in a freshwater tidal march.
The American Midland Naturalist, 143(1), 53-63.
Chefaoui, R.M., Hortal, J. & Lobo, J.M. (2005). Potential distribution modelling, niche characterization and conservation status assessment using GIS tools: a case study of Iberian Copris species. Biological
Conservation, 122, 2, 327-338.
Doude van Troostwijk, W. J. (1978). Muskrat control in the Netherlands. Wildlife Division, Ministry of
Agriculture and Fisheries.
Ecke, F., Henry, A. & Danell, K. (2013). Landscape-based prediction of the occurrence of the invasive muskrat (Ondatra zibethicus). Finnish zoological and Botanical Publishing Board, 51(3), 325-334.
Galparsoro, I., Borja, A., Bald, J., Liria, P. & Chust, G. (2009). Predicting suitable habitat for the European lobster (Homarus Gammarus), on the Basque continental shelf (Bay of Biscay), using Ecological-Niche Factor Analysis. Ecological Modelling, 220, 4, 556-567.
Heidinga, D. (2006). Pluizige plaagdieren, Ecologie en bestrijding van de muskusrat. Wetenschapswinkel Biologie, rapport 70.
Hirzel, A. H., Hausser, J., Chessel, D. & Perrin N. (2002). Ecological-niche factor analysis: how to compute habitat suitability maps without absence data? Ecology, 83(7), 2027-2036.
LCCM (2007). Landelijk jaarverslag 2006 muskus- en beverrattenbestrijding.
Unie van Waterschappen (2016). Landelijk Jaarverslag Muskus- en Beverratten 2015.
Van Loon, E. E., Bos, D. Van Hellenberg Hubar, C. J. & Ydenberg R. C. (2016). A historical perspective on the effects of trapping and controlling the muskrat (Ondatra zibethicus) in the Netherlands.
Skyrienne, G. & Paulauskas, A. (2012). Distribution of invasive muskrats (Ondatra zibethicus). Ekologija, 58(3), 357-367.
Xuezhi, W., Weihua, X., Zhiyun, O., Jianguo, L., Yi, X., Youping, C., Lianjun, Z. & Junzhong, H. (2008). Application of ecological-niche factor analysis in habitat assessment of giant pandas. Acta Ecologica Sinica, 28(2), 821-828.
Appendix 1. GIS-commands
Data: Gemeentegrenzen 2013
Source: ArcGIS online
Creating DEM of gemeente Schermer by using clickable maps
Importing maps
Select By Attributes>”GM_NAAM” = ‘Schermer’>Gemeente grezen>OK>Selection>Create layer from selected features.
Data: AHN2 (clickable map). Source: geoportal of university of Amsterdam:
http://geodata.science.uva.nl/UvAGeodata/ClickableMaps/World/Europe/Netherlands/AHNClicka bleMap/AHNClickableProvinces.html
DEM of The Netherlands, used data:
• I19gn1.tif • i19ez1.tif • I19dn1.tif • I19dn2.tif • I19bz2.tif • I19bz1.tif • i19dn2.tif
Mosaic layers into one new raster
Search>Mosaic to new raster>Input Features: i19gn1.tif, i19ez1.tif, i19dn1.tif, i19dn2.tif, i19bz2.tif, i19bz1.tif, i19dn2.tif >PixelType:32_BIT_FLOAT>number of bands:1>OK
Creating aspect and slope map
Aspect map of DEM-layer “Schermer_DEM”
Search>aspect>input: “Schermer”>OK
Slope map of DEM “Schermer_DEM”
Search>slope>input:”Schermer”>OK
Source: Waterschap Hoogheemraadschap Hollands Noorderkwartier
Used data:
BegroeidTerreindeel_v Waterdeel_v
Obtaining research area of 25m around water bodies
Clipping ‘waterdeel_v’ into ‘Schermer_waterdeel’
Geoprocessing>Clip>Input Feature: “Schermer”> Clip features: ‘waterdeel_v’>OK
Selecting research area of 25m around water bodies. Output: ‘Schermer_researcharea’
Search>Buffer>Input Features: ‘Schermer_waterdeel’>distance: linear unit: 25m> output dissolve type: all>OK
Extract mask layer of aspect and ‘Schermer_researcharea’
Search> extract by mask> Input raster: ‘Schermer_aspect’>Input raster or feature mask data: ‘Schermer_researcharea’>OK
Change cell size so that maps can be imported to R-studio
Search>resample>’aspect_mask’>X:5m>Y:5m>OK Creating rasters from polygon data
Changing polygon to raster ‘BegroeidTerreindeel_V’
Search>polygon to raster> input features: ‘BegroeidTerreindeel_V’>Value fied:
fysiekVoorkomen> cellsize:5 (arcmap is not able to perform raster with smaller cellsize)
Data: DNA.gdb
Source: Waterschap Hoogheemraadschap Hollands Noorderkwartier. Dataset was developed for research about possibility of using DNA tests in water to increase muskrat catchments.
Used data: Grondsoorten
Changing polygon to raster ‘Grondsoorten’
Search>polygon to raster> input features: ‘Grondsoorten’>Value fied: fysiekVoorkomen>cellsize: 5
Soil map of 5mx5m of research area of 20m around water bodies
Search> extract by mask> Input raster: ‘raster_soil’>Input raster or feature mask data: ‘Aspect_mask_resample’>OK
Terrain feature map of 5mx5m of research area of 20m around water bodies
Search> extract by mask> Input raster: ‘raster_bt’>Input raster or feature mask data: ‘Aspect_mask_resample’>OK
Slope map of 5mx5m of research area of 20m around water bodies
Search> extract by mask> Input raster: ‘Schermer_slope’>Input raster or feature mask data: ‘Aspect_mask_resample’>OK
Export data
begroeid_terrein_mask_5m>Data>Export Data>NoData as: NaN>Save soil_mask_5m>Data>Export Data>NoData as: NaN>Save
aspect_mask_5m>Data>Export Data>NoData as: NaN>Save slope_mask_5m>Data>Export Data>NoData as: NaN>Save
Point data to binary-map
Selecting from attribute table Selection>Select By Attributes>
MUSKUSRAT_RAM_OUD >0 OR MUSKUSRAT_MOER_OUD>0 OR MUSKUSRAT_RAM_JONG>0 OR MUSKUSRAT_MOER_JONG>0 AND SOORT_COMBINATIE = 1 (selectie conibear trap)
>OK>VangstRegistratie_NoordHolland>Selection>Create layer from selected features
Selection “Schermer”
Geoprocessing>Clip>input:”muskrat”> Clip features: “Schermer’’>OK Creating boolean map
‘muskrat_schermer_coni’>open attribute table>options>add field>name: ‘booleannumber’>ok>calculator field booleannumber>’booleannumber’=1.
‘Schermer’>open attribute table>options>add field…>name>booleannumber>OK>field calculator>booleannumber:0>OK
Search>polygon to raster>input:’Schermer’>value field: booleannumber>Cellsize: 5m>OK
Mosaic To New Raster> Input Rasters: ‘muskrat_schermer_coni’ ‘Schermer_polygontoraster’> ellsize: catchment_point to raster>pixel type: 16_bit_unsigned>mosaic operator:maximum>OK
‘binarymap’>properties>unique values>apply
Binary map of cell size: 5mx5m of research area of 20m around water bodies
Search> extract by mask> Input raster: ‘Schermer_slope’>Input raster or feature mask data: ‘Aspect_mask_resample’>OK
Appendix 4. Script RStudio
# Ecological-niche factor analysis on muskrats in Schermer, Noord Holland # Bachelor thesis Thirza van Noppen
#Install packages install.packages('raster') install.packages('rgdal') install.packages('adehabitatHS') install.packages("randtests") library(raster) library(ade4) library(adehabitatHS) library(randtests) setwd("G:/area25m")
#Importing "vangstdata" muskrats, change extent and convert into single vector pres <- raster("Boolean.tif")$Boolean
extent(pres)<-extent(c(111398,123798,507801,517801)) presence<-as.data.frame(pres)
#change NA to 0
presence[is.na(presence)]<- 0 #Importing raster maps into R aspect <-raster("Aspect.tif") slope <-raster("Slope.tif") soil_1 <-raster("Soil.tif") soil_2 <-raster("Soil.tif") soil_3 <-raster("Soil.tif") soil_4 <-raster("Soil.tif") beg_ter_1 <-raster("Begr_ter.tif") beg_ter_2 <-raster("Begr_ter.tif") #Change extent of raster layers
extent(soil_1)<-extent(c(111398,123798,507801,517801)) extent(soil_2)<-extent(c(111398,123798,507801,517801)) extent(soil_3)<-extent(c(111398,123798,507801,517801)) extent(soil_4)<-extent(c(111398,123798,507801,517801)) extent(beg_ter_1)<-extent(c(111398,123798,507801,517801)) extent(beg_ter_2)<-extent(c(111398,123798,507801,517801)) #Find clay soils and replace by 1, leave rest 0
clay<-soil_1 clay[clay==7]<-100
clay[clay==10]<-100 clay[clay<11]<-0 clay[clay>11]<-1
#Find bog soils and replace by 1, leave rest 0 bog<-soil_2
bog[bog==8]<-100 bog[bog<11]<-0 bog[bog>11]<-1
#Find sand soils and replace by 1, leave rest 0 sand<-soil_3
sand[sand==5]<-100 sand[sand<11]<-0 sand[sand>11]<-1
#Find zavel soils and replace by 1, leave rest 0 zavel<-soil_4
zavel[zavel==1]<-100 zavel[zavel==2]<-100 zavel[zavel<11]<-0 zavel[zavel>11]<-1
#Find trees(shadow) and replace by 1, leave rest 0 trees<-beg_ter_1
trees[trees==8]<-100 trees[trees==29]<-100 trees[trees<30]<-0 trees[trees>30]<-1
#Find reed and replace by 1, leave rest 0 reed<-beg_ter_2 reed[reed==18]<-100 reed[reed<30]<-0 reed[reed>30]<-1 var<-stack(aspect,slope,clay,bog,sand,zavel,trees,reed) spdf<-as(var, "SpatialPixelsDataFrame")
#Convert data into dataframe var1 <- as.data.frame(aspect) var2 <- as.data.frame(slope) var3 <- as.data.frame(clay) var4 <- as.data.frame(bog)
var8 <- as.data.frame(reed)
map_nan = cbind(var1, var2, var3, var4,var5, var6,var7,var8,presence) #Remove NaN, transform back to list
map <- map_nan[complete.cases(map_nan), ] map_list <-as.list(map[1:8])
# create spatial points data frame spg <- map_list
# Species map species <- map[,9]
#Visualization data, kernel density #Aspect muskrat <-which(species==1) d_Aspect<-density(map_list$Aspect) d_m_Aspect<-density(map_list$Aspect[muskrat]) par(mar = c(5, 5, 3, 5)) plot(d_Aspect,xlim = c(0,400),ylim = c(0,0.006),col=c("red"),main="",xlab="",ylab="Density",yaxs="i",xaxs="i") par(new=TRUE)
plot(d_m_Aspect, ylab="",xlim = c(0,400),ylim=c(0,0.006),col=c("blue"),main="Kernel density plot of aspect",xlab="Degrees",xaxt = "n", yaxt = "n",yaxs="i",xaxs="i",lty = 2)
legend("topright",lwd=1,c("Total research area, B=5.971","Nest-sites, B=55.55"),col=c("red","blue")) #Slope d_Slope<-density(map_list$Slope) d_m_Slope<-density(map_list$Slope[muskrat]) par(mar = c(5, 5, 3, 5)) plot(d_Slope,xlim = c(0,70),ylim = c(0,0.45),col=c("red"),main="",xlab="",ylab="Density",yaxs="i",xaxs="i") par(new=TRUE)
plot(d_m_Slope, ylab="",xlim = c(0,70),ylim=c(0,0.45),col=c("blue"),main="Kernel density plot of slope",xlab="Degrees",xaxt = "n", yaxt = "n",yaxs="i",xaxs="i",lty = 2)
legend("topright",lwd=1,c("Total research area, B=0.07926","Nest-sites, B=0.6138"),col=c("red","blue")) #Clay
d_clay<-density(map_list$Soil)
d_m_clay<-density(map_list$Soil[muskrat]) par(mar = c(5, 5, 3, 5))
plot(d_clay,xlim = c(0,1),ylim = c(0,10),col=c("red"),main="",xlab="",ylab="Density total research area",xaxt = "n",yaxs="i",xaxs="i")
plot(d_m_clay, xlim = c(0,1), ylab="",ylim=c(0,2),col=c("blue"),main="Kernel density plot of clay soil",xlab="Percentage clay soil",xaxt = "n", yaxt = "n",yaxs="i",xaxs="i",lty = 2)
axis(1,at=seq(0, 1, by=0.2),labels=c("0","20","40","60","80","100")) axis(side=4)
mtext("Density nest-sites",side=4,line=3)
legend("topright",lwd=1,c("Total research area, B=0.02805","Nest-sites, B=0.2436"),col=c("red","blue")) #Bog
d_bog<-density(map_list$Soil.1)
d_m_bog<-density(map_list$Soil.1[muskrat]) par(mar = c(5, 5, 3, 5))
plot(d_bog,xlim = c(0,1),ylim = c(0,13),col=c("red"),main="",xlab="",ylab="Density total research area",xaxt = "n",yaxs="i",xaxs="i")
par(new=TRUE)
plot(d_m_bog, xlim = c(0,1), ylab="",ylim=c(0,4),col=c("blue"),main="Kernel density plot of bog soil",xlab="Percentage bog soil",xaxt = "n", yaxt = "n",yaxs="i",xaxs="i",lty = 2)
axis(1,at=seq(0, 1, by=0.2),labels=c("0","20","40","60","80","100")) axis(side=4)
mtext("Density nest-sites",side=4,line=3)
legend("topright",lwd=1,c("Total research area, B=0.02426","Nest-sites, B=0.1068"),col=c("red","blue")) #sand
d_sand<-density(map_list$Soil.2)
d_m_sand<-density(map_list$Soil.2[muskrat]) par(mar = c(5, 5, 3, 5))
plot(d_sand,xlim = c(0,1),ylim = c(0,70),col=c("red"),main="",xlab="",ylab="Density total research area",xaxt = "n",yaxs="i",xaxs="i")
par(new=TRUE)
plot(d_m_sand, xlim = c(0,1), ylab="",ylim=c(0,1),col=c("blue"),main="Kernel density plot of sandy soil",xlab="Percentage sandy soil",xaxt = "n", yaxt = "n",yaxs="i",xaxs="i",lty = 2)
axis(1,at=seq(0, 1, by=0.2),labels=c("0","20","40","60","80","100")) axis(side=4)
mtext("Density nest-sites",side=4,line=3)
legend("topright",lwd=1,c("Total research area, B=0.005521","Nest-sites, B=0.4896"),col=c("red","blue")) #Zavel
d_zavel<-density(map_list$Soil.3)
d_m_zavel<-density(map_list$Soil.3[muskrat]) par(mar = c(5, 5, 3, 5))
plot(d_zavel,xlim = c(0,1),ylim = c(0,14),col=c("red"),main="",xlab="",ylab="Density total research area",xaxt = "n",yaxs="i",xaxs="i")
par(new=TRUE)
plot(d_m_zavel, xlim = c(0,1), ylab="",ylim=c(0,1.4),col=c("blue"),main="Kernel density plot of fine sand soil",xlab="Percentage fine sand soil",xaxt = "n", yaxt = "n",yaxs="i",xaxs="i",lty = 2)
#Trees
d_t<-density(map_list$Begr_ter)
d_m_t<-density(map_list$Begr_ter[muskrat]) par(mar = c(5, 5, 3, 5))
plot(d_t,xlim = c(0,1),ylim = c(0,150),col=c("red"),main="",xlab="",ylab="Density total research area",xaxt = "n",yaxs="i",xaxs="i")
par(new=TRUE)
plot(d_m_t, xlim = c(0,1), ylab="",ylim=c(0,1),col=c("blue"),main="Kernel density plot of trees",xlab="Percentage trees",xaxt = "n", yaxt = "n",yaxs="i",xaxs="i",,lty = 2)
axis(1,at=seq(0, 1, by=0.2),labels=c("0","20","40","60","80","100")) axis(side=4)
mtext("Density nest-sites",side=4,line=3)
legend("topright",lwd=1,c("Total research area, B=0.002572","Nest-sites, B=0.2482"),col=c("red","blue")) #Reed
d_r<-density(map_list$Begr_ter.1)
d_m_r<-density(map_list$Begr_ter.1[muskrat]) par(mar = c(5, 5, 3, 5))
plot(d_r,xlim = c(0,1),ylim = c(0,55),col=c("red"),main="",xlab="",ylab="Density total research area",xaxt = "n",yaxs="i",xaxs="i")
par(new=TRUE)
plot(d_m_r, xlim = c(0,1), ylab="",ylim=c(0,4),col=c("blue"),main="Kernel density plot of reed",xlab="Percentage reed",xaxt = "n", yaxt = "n",yaxs="i",xaxs="i",lty = 2)
axis(1,at=seq(0, 1, by=0.2),labels=c("0","20","40","60","80","100")) axis(side=4)
mtext("Density nest-sites",side=4,line=3)
legend("topright",lwd=1,c("Total research area, B=0.007263","Nest-sites, B=0.1068"),col=c("red","blue"))
#pca
pc <- dudi.pca(map_list,scannf=FALSE,nf = 8) #ENFA
# enfa_output <- enfa(pc, presence,scannf=FALSE,nf = 4) enfa_output <- enfa(pc, species,scannf=FALSE,nf = 8) scatter(pc)
hist(enfa_output)
hist(enfa_output, scores = FALSE, type = "l")
## scatterplot scatter(enfa_output)
#renfa <- randtest(enfa_output) #plot(renfa)
# Produce habitat suitability map # install package sp
library(sp)
