Landscape ecology and wild rabbit (Oryctolagus cuniculus) habitat modeling in the Mediterranean region

in 277–283. Abstract Landscape ecology and wild rabbit ( Oryctolagus cuniculus ) habitat modeling in the Mediterranean region. — Landscape modification is one of the reasons for the decrease in rabbit populations. The objective of this study was to model wild rabbit habitat using landscape ecology to create a diagnosis method able to assess habitat quality at a large scale. Rabbit presence/absence was recorded on 536 plots of 1 ha. Spotlight transect counts indicated a low relative abundance (KIA = 2.3 rabbits/km). We produced a land use map with metric precision using remote sensing. Water, bare soil, herbaceous, shrubs and trees were identified. Landscape structure and diversity were evaluated using variables available in FRAGSTATS. A logistic regression was performed to assess the link between rabbit presence/absence and landscape structure. Our results indicate that a suitable habitat has a high diversity, a medium number of patches and a small proportion of shrubs. These results could be used to diagnose the landscape prior to any management action to enhance rabbit populations and conversely be helpful as a tool of integrated control in the cases of local outbreaks with agricultural damages.


Introduction
Populations of European wild rabbit (Oryctolagus cuniculus) have decreased in France since the 1950s (Marchandeau, 2000).This decline is mainly due to diseases (myxomatosis and rabbit viral haemorrhagic disease), evolution of agricultural practices and excessive hunting (Trout et al., 1992;Villafuerte et al., 1995;Marchandeau et al., 1998).Rabbits are considered a keystone species in the Mediterranean region (Delibes-Mateoset al., 2007).Because they are prey for many threatened species such as Bonelli's eagle Aquila fasciata or eagle owl Bubo bubo (Delibes & Hiraldo, 1981) recovery of rabbit populations could have a positive effect on the conservation of these species (Moreno & Villafuerte, 1995).Furthermore, the rabbit was the principal game species in France until the 1950s and it continues to be a main target species To enhance rabbit populations, detailed knowledge of habitat requirements is fundamental.In recent years, many studies have been carried out along these lines (Rogers & Myers, 1979;Calvete et al., 2004;Carvalho & Gomes, 2004;Beja et al., 2007;Serrano Pérez et al., 2008) and several reports have pointed out the importance of the landscape structure.More specifically, they have observed that the interspersion of shelter and open areas is linked to the predation risk (Palomares & Delibes, 1997;Villafuerte & Moreno, 1997).
The aim of this study was to create a diagnostic method able to assess habitat quality at a finer spatial scale than in previous works and in large areas (several hundred or thousand hectares).This could allow us to develop guidelines for managers and help them to increase rabbit populations.We first developed a high definition (metric precision) land-use map using remote sensing and automated photo-interpretation (Lucas et al., 2007).We then analysed the landscape structure to create a habitat suitability model in order to identify factors related to rabbit distribution.

Study area
The study area is located in the southeast of France and covers the Natura 2000 (ZPS) area of Durance valley (figs.1A, 1B).It is 240 km long and 2 km wide (1 km on both sides of the river) for an area of 46,500 ha.The altitude varies from 15 m to 780 m a.s.l. and there is a Mediterranean climate (average annual temperature and rainfall, 12.4°C and 720 mm respectively) (Benichou & Le Breton, 1987).The land use in the study area is: 12% of artificial area, 38% of agricultural area, 42% of forest and semi-natural area and 8% of wetland and water (Union européenne-SOeS, 2006).The main natural areas are riparian forest with Populus sp., Alnus sp., Salix sp. and Quercus sp.Agricultural practices are mainly mixed farming with crop, arboriculture and breeding.
Six pilot areas of about 100 ha each were demarcated and divided into plots of 1 ha.Using spotlight transect counts, an average kilometric abundance index of 2.3 rabbits/km was obtained (Ballinger & Morgan, 2002).It indicated a low relative abundance of rabbit.

Determination of presence/absence areas
Rabbit presence/absence was determined on 536 plots of 1 ha distributed in the various pilot areas.The plot area was close to the size of the home range of the rabbit (Lombardi et al., 2007;Marchandeau et al., 2007).Each plot was randomly prospected during 15' noting each sign of rabbit presence: pellets, latrines, burrows and scratches.Two field sessions were conducted during the third quarter of 2008 and 2009 and two others during the first quarter of 2009 and 2010.Thus, each plot was prospected four times.The information was synthesized by considering that rabbits were present on a plot if their presence was noted at least once.Rabbits were considered absent if their absence was noted four times on a plot.

Land use cartography
In order to quantify landscape structure by calculating several metrics, a land use map with metric precision was performed using remote sensing technology.Aerial photographs with dot pitch of 50 cm were taken in July 2009 using the near infrared channel and the red, green and blue channel.An automated photointerpretation was made in two steps with eCognition (Trimble) software (Baatz et al., 2001) using the object-oriented classification method which is suitable for very high resolution (Xiaoxia et al., 2005).First, we performed a segmentation consisting of delimiting polygons with homogenous land use.Each polygon was then classified as one of the land use types.The five different land use types identified were: water, bare soil, herbaceous (< 1 m), shrubs (1 << 7 m) and trees (> 7 m).Finally, the data were rasterized with cells of 4 m² (2 m x 2 m) (Xiaoxia et al., 2005).

Calculation of landscape metrics
Landscape structure and diversity were evaluated using 20 variables.Eleven were calculated using ArcMap™ 9.1 software using FRAGSTATS interface (McGarigal & Marks, 1994)

Fig. 1. Área de estudio en el sureste de Francia (A) y localización de las seis áreas experimentales (B).
randomly positioned in cover areas and the nearest open area.For each plot of 1 ha, the surface of cover was calculated and one point was used for 100 m² of cover with a minimum of 30 points.The same was applied for the NEARESTCOVER which indicates the average mean distance to the nearest cover.All the variables are quantitative except PR which is qualitative with three modalities (one or two patch types: PR = a, three patch types: PR = b, four or five patch types: PR = c).

Statistical analysis
All the statistical analysis were performed using R© software (2.12.2 version) (R Development Core Team, 2011).The descriptive analysis reveals the presence of plots with very high values of NUMP and AWMSI.The specific study of these plots indicates that they were composed of orchards and row crops.Because of this particular landscape structure, remote sensing was poor.Thus, they were considered as outliers and removed from the data set which was finally composed of 516 measures.The correlation matrix indicates some colinearity and this was taken into account by computing the Variance Inflation Factor (VIF) with a threshold of 3 to test colinearity among the variables included in the different models.
First, we compared the means of the variables from presence versus absence areas.Homoscedasticity was tested by an F-test.In cases of homoscedasticity, a t-test was performed.Otherwise, a Welch's test was carried out.
A generalized linear model was then used with Logit as a link function to assess the link between rabbit presence/absence (binary independent variable) and the landscape metrics (dependant variables) (Traissac et al., 1999).The data set was randomly divided in two parts: 80% of the sample for model calibration and the remaining 20% for validation.
The first selection was done by minimizing Akaike Information Criterion (AIC) (Akaike, 1974;Lebreton et al., 1992;Posada & Buckley, 2004) with a stepwise selection starting with null model (function 'step' under R software with argument direction = 'both') (R Development Core Team, 2011).It consists of adding variables one by one and trying to remove those of the previous steps in order to minimize AIC.On the first model obtained by stepwise selection, Wald's test was used to identify variables whose coefficient was not significantly different from zero (Hauck & Donner, 1977); these variables were removed.Finally, the sign of the coefficient of the different variables was examined.When we failed to find a logical and biological explanation (inconsistencies with other results in this study and with literature, no colinearity), the variable was removed.Interactions of first degree were tested with Wald´s test in the model obtained by stepwise selection and in the final model.The goodness of fit was tested with Hosmer-Lemeshow test (Lemeshow & Hosmer, 1982;Hosmer & Lemeshow, 2000) and McFadden's pseudo r-squared.The discrimination capacity of the model was evaluated by means of a confusion matrix and Area Under the Receiver Operating Characteristic (ROC) Curve (Agresti, 2002).

Results
Mean difference between presence and absence areas are summarized in table 1.A significant difference was found for 15 variables.Seven variables had a higher value in the presence areas: IJI, NUMP, ED, AWMSI, AWMPFD, PR and PROPTREES.In contrast, MNN, MPS, TCA, CAD, LPI, PROPSHRUBS, NEARESTCO-VER and NEARESTOPENAREA were significantly higher in areas of absence.
The first model obtained by stepwise selection was composed of the following variables: NUMP, NUMP², PR, PROPSHRUBS, IJI and TCA.The VIF indicates that there was no colinearity among these variables.According to the Wald test, the coefficient of IJI was not significantly different from 0 (P = 0.151), so this variable was removed from the model.The positive coefficient of TCA seems to specify that large patches with an important core area would be more suitable for rabbit.This result was not due to colinearity and it was not consistent with other results of the present study or with the literature; so this variable was also removed.The final model was made up of four variables: NUMP,  2).All of them were significant according to the Wald's test.Interactions were tested in the various previous models but were not significant.Figure 2 (A and B) illustrates the effect of the variables on the probability of rabbit presence.Figure 2A indicates that the probability of presence is higher when there are three or more patch types in the landscape, suggesting it should be diversified to be favourable.Figure 2B summarises the effect of the proportion of shrubs and the number of patches on the probability of presence.The probability decreased when the proportion of shrubs increased.Because of the presence of the square term of NUMP, the probability of presence correlated positively with the number of patches up to about 100 patches/ha, after which it was negatively correlated.It thus reaches its maximum when there are about 100 patches per hectare.The Hosmer-Lemeshow test (3.99,P = 0.86) indicated that the model correctly fitted the data.McFadden's pseudo r-square (0.19) was also computed and indicated an adequate fit.The Variance Inflation Factor revealed no colinearity.Using the validation data set, the area under the ROC curve (AUC) was 0.79.This value shows that the model had a high capacity for discrimination.Sensitivity and the specificity were calculated using the confusion matrix and were 76% and 73%, respectively.

Discussion
The aim of this study was to identify landscape factors that influence rabbit distribution.Climatic conditions and altitude were thus not used as predictor variables.However, even if there is a large range of altitude on the study area, climatic conditions are consistent with rabbit requirements.To detect the presence of rabbits, all the signs of presence were taken into account rather than only the traditionally used pellets.The risk of declaring an absence where the rabbit was present was therefore minimal.The endpoint was to identify suitable habitats on the basis of rabbit distribution.Their distribution also depends on other factors, however, such as disease impact, and it can vary from one year to the next.To minimize these biases, four field sessions were conducted, spread over three years and an area was considered as suitable if rabbit presence was recorded at least once.However, the absence of rabbit from some suitable areas cannot be totally ruled out.
The probability threshold used to discriminate rabbit presence and absence was set at 0.75.This value permitted to maximise the discrimination capacity of the model.It only discriminated between suitable and unsuitable habitats because it works as a binary variable.However, when the model is used for prediction, it can be modified according to the risk to be minimised.Increasing the threshold would limit the number of false positive cases and decreasing it would limit false negative cases.
In this study, the land use map was created using remote sensing which combines high definition aerial photography with automated photo-interpretation. Aerial photography with a dot pitch of 50 cm allowed us to make a more precise map than those using classical aerial photography or 1:25,000 maps.Thus, the map was composed of pixels of 4 m² which seem to be consistent with a rabbit habitat study.This high precision could be interesting to precisely evaluate habitat structure in relation to the perception scale of the rabbit.Moreover, the use of remote sensing can work on very large areas.This technique avoids manual digitalisation which is very time consuming on areas of several hundred hectares.However, it has a disadvantage because it only takes into account the highest vegetation layer.For example, it ignores the presence of herbaceous or shrubs under trees.
The results presented in table 1 suggest that fifteen variables could be useful to describe rabbit habitat,  The variables IJI, AWMSI and AWMPFD underline the importance of the fractal dimension (Jimenez et al., 2006) and interspersion (Carvalho & Gomes, 2004) between patches.
The statistical significance of the variables NUMP, PROPSHRUBS and NEARESTCOVER indicates that a suitable habitat is a patchwork with some evenly distributed shrub cover.This is consistent with the results of Moreno & Villafuerte (1995) that indicated that rabbit abundance was greatest at less than 20 m from scrub cover.The importance of the arrangement of cover and open areas can be considered an antipredator strategy.Cover is safe for the rabbit to hide from the birds of prey during the day.On the contrary, open areas are safer at night because carnivorous mammals use cover for hunting (Moreno et al., 1996).
This study provides some characteristics of rabbit habitat obtained using remote sensing using a fine spatial scale.Based on our results, we have a diagnostic method which could be used to assess habitat quality concerning rabbit requirements.This process could enable accurate discrimination between suitable and unsuitable areas over areas of several hundred or thousand hectares.Evaluation of unsuitable areas might identify variables of inadequate value.Such information would be of valuable help for managers trying to improve habitat quality in order to enhance rabbit populations.Conversely, this method could be used as a tool of integrated control to reduce rabbit populations (Boag, 1987).In the case of local outbreaks causing agricultural damage, modification of habitat structure could create an unsuitable habitat and lead to a long-term decrease of rabbit populations.
Fig. 1.Study area in the southeast of France (A) and location of the six pilot areas (B).

Fig. 2 .
Fig. 2. Relationship between probability of presence and patch richness (A) and relationship between probability of presence, proportion of shrubs and number of patches for PR = b (B).

Table 1 .
Mean difference between variables in absence areas (Aa) vs. presence areas (Pa).(Forabbreviations of variables see Material and methods.Tabla 1. Diferencias entre las medias de las variables en las zonas de ausencia (Aa) vs. presencia (Pa).(paralas abreviaturasde las variables, ver Material and methods.)

Table 2 .
Coefficients of the variables retained in the model and P-value of Wald test.