Introduction

The white-tailed deer Odocoileus virginianus is the most important game species in Mexico (^{Ortega-Santos et al. 2011}). This deer species can be found in different habitats including dry tropical forests, temperate forests, and xeric shrubland, among others (^{Gallina et al. 2010}). In addition, this species persists in farming and ranching areas interspersed with patches of the original habitat (^{Mandujano et al. 2014}). However, white-tailed deer populations are affected by poaching, irrational use, habitat disruption, competition for resources with domestic livestock, and incidence of parasites (^{Ortega-Santos et al. 2011}). A proper utilization of this species requires estimating population density and its relationship with different features of the habitat, vegetation and landscape, among others (^{Ortiz-Martinez et al. 2005}; ^{López-Tellez et al. 2007}; ^{Medina-Torres et al. 2008}; ^{Mandujano et al. 2013}; ^{Ramos-Robles et al. 2013}; ^{Barrera-Salazar et al. 2015}; ^{García-Marmolejo et al. 2015}). Unlike other regions of the country — for example the northeast (^{Villarreal 1999}), La Michilía in Durango (^{Galindo-Leal and Weber 1998}), the Mixteca Poblana (^{Villarreal-Espino 2006}), and various regions of the state of Oaxaca (^{Mandujano 2016}) — management of the white-tailed deer is not an extensive practice in the state of Veracruz; hence, a larger number of studies are required for its populations (^{Mandujano 2011}). Investigations conducted on this species in Veracruz include topics such as the distribution of deer, habitat quality and management in Management Units for Wildlife Conservation and Sustainable Use (UMAs) (^{Bello and Mandujano 1994}; ^{Gallina et al. 2007}; ^{Delfín-Alfonso et al. 2009}; ^{Bolívar-Cime and Gallina 2013}, ^{Serna-Lagunes et al. 2013a}, ^b; ^{Contreras-Verteramo et al. 2014}; ^{Mandujano and López-Tello 2015}).

The objective of this work was to estimate and compare the population density of the white-tailed deer in five locations and its relationship with variables in the different types of vegetation in a landscape of the municipality of Tantoyuca, northern Veracruz. Ethnic groups such as the Huasteco and the Nahuatl coexist in this mountainous area of the Huasteca Alta Veracruzana (^{Reyes 2009}); 77 % of these human populations have been classified as living in poverty, under marginality conditions and under a high social gap (^{SEFIPLAN 2016}). The main productive activities in the region are temporal agriculture and livestock raising. Similar to many regions of the state, the white-tailed deer is subjected to uncontrolled hunting. In some cases, local inhabitants complain that the deer damages crops, although there are also communities that protect and hunt it for their own consumption. However, to date there are no UMAs to exploit it legally. Consequently, the present work also aims to issue some recommendations for the protection and potential management of the white-tailed deer and its habitat in this region.

Materials and Methods

Study area. The study was conducted in five locations (Chiquero, Pensador Mexicano, Porvenir, Monte Grande and Corralillo) located in the municipality of Tantoyuca in the northern portion of the state of Veracruz (Figure 1a). The municipality is located at 21° 06’ to 21° 40’ N and -97° 59’ to -98° 24’ E at an altitudinal range of 10 to 300 m. Climate is warm subhumid with a mean annual precipitation of 1200 mm and a mean annual temperature of 24 °C (^{INEGI 2009}). The landscape lies within the Calabozo river basin, a tributary of Panuco river, and in small streams emptying into Topila estuary (^{INEGI 2009}). The municipality has an area of 1,303 km² and comprises 511 rural and one urban localities. The original representative vegetation in this region was sub-evergreen forest with patches of high tropical forest (^{Rzedowski 1978}). Current land use in the municipality includes agriculture (13 %), pastureland (57 %), secondary vegetation (28 %) and urban area (< 1 %). Local economic activities include agriculture, livestock raising (cattle, pigs and sheep) and poultry (^{INEGI 2009}). Consequently, the landscape is a mosaic of different vegetation types with varying degrees of disturbance, alternating with small patches of the original forest, pastureland and cultivated areas (Figure 1b).

Figure 1 (a) Geographic location and dominant vegetation types according to ^{INEGI (20005)} in the region of study, municipality of Tantoyuca, northern state of Veracruz. (b) Landscape heterogeneity. 

Field sampling. The five study localities were selected based on ease of field sampling and the interest of the local community in caring for and making good use of deer within their premises. Vegetation was characterized using the same protocol as in previous studies (^{Mandujano et al. 2013}; ^{Ramos-Robles et al. 2013}; ^{Barrera-Salazar et al. 2015}). Basically, eight 500-m transects were set in each location, with 11 sampling points separated by 50 m between them in each location, for a total of 440 points. Vegetation was sampled using the point quadrat method, also used in previous studies on deer. The variables studied were composition of tree species, richness (N0), Shannon diversity (H), Hill’s diversity (N1), evenness (E1), density (Arb_den), height (Arb_alt) and basal area (Arb_ab) of the tree stratum; understory species richness (S); and protection cover at 0-50 cm (cob_50) and 51-100 cm (cob_100). N0, H, N1 and E1 were calculated using the package R (^{Bocard et al. 2010}). Samplings were carried out between April and June 2014.

For deer pellet counts, 40 strip transects measuring 500 x 2 m sere set, with a minimum separation of 500 m between them. This distance attempted to prevent that the same animal leave traces in different transects; however, this might not be true in all cases (^{Mandujano 2014a}). Each transect was divided into 10 sampling plots measuring 50 x 2 m for counting and collecting all deer pellet groups, taking care not to miss feces within the transect, nor include those observed outside the strip, to avoid any bias. Samplings were conducted during the dry season from April to June 2014. All transects were sampled by the same person assisted by local inhabitants, in order to reduce any observational bias during sampling. Location data were collected using GPS, and were subsequently imported to the program QGIS 2.4 (http://qgis.org/es/site/) for its spatial projection.

Classification of vegetation. Each transect was classified according to the dominant vegetation type in INEGI’s National Forest Inventory map (2005) (Figure 1a). However, considering the heterogeneity of the landscape, a single transect frequently covered one or more vegetation types; in these cases, transects were classified according to the dominant vegetation. In addition to the classification of ^{INEGI (2005)}, in this study we classified the 40 transects according to their similarity in the composition of tree species. To this end, hierarchical clustering methods were used for species abundance data through multivariate clustering analysis. The original matrix was 40 transects by 101 species. However, in order to facilitate the data analysis and interpretation, 57 species with less than 15 individuals were removed. The analysis used the procedures described in ^{Bocard et al. (2010)} and the ade4 gclus, vegan, and cluster packages in R. The matrix of Euclidean distances between taxonomic groups was calculated using the hclust function. The classification methods used were the unweighted-pair-group method of analysis (UPGMA) and Ward’s minimum variance. The results of hierarchical clustering were interpreted and compared using some indices and plots provided by these packages. Finally, the species assigned to groups or “clusters” were examined with the vegemite function of the vegan package (^{Bocard et al. 2010}).

Once the classification of transects in each group was completed, the vegetation variables (N0, H, E1, Arb_den, Arb_alt, Arb_ab, S, cob_50 and cob_100) were tested for statistical differences using one-way ANOVAs. Previously, the Shapiro test was run to confirm the normality of the residuals, and Bartlett’s test to determine the homogeneity of variances. In case of significance in the analysis of variance, Tukey’s HSD test was used a posteriori for determining differences between groups (^{Crawley 2013}).

Relationship between density and vegetation. Population density was estimated with the program PELLET version 2.0 for Excel (^{Mandujano 2014b}). The equation for estimating population density was: D= (Dpg*Np)/(Tdf*Tdep), where: D = deer population density (number of individuals per km²), D_pg = average density of pellet groups per transect, N_p = Number of transects sampled, T_df = defecation rate, and T_dep = days of deposition. This procedure uses three defecation rates (13, 19 and 26 pellet groups/ind/day) obtained in several studies, while exposure rate is calculated from the latest rainfall in the study site and the sampling date in each transect. Further details on this procedure are provided in the work of ^{Mandujano (2014a}, ^b).

Prior to the statistical comparison of deer density, preliminary analyses of the data were carried out using the quantil-quantil plot function in R, where theoretical and sampling residuals were plotted (^{Crawley 2013}). The influence function was used to detect any observations (“outliers”) that might skew to a greater extent the density estimates. Deer density was estimated at four levels: 1) for the whole municipality, by grouping all transects; 2) for each of the five locations; 3) by vegetation type according to the INEGI classification (^{INEGI 2005}); and 4) by groups or clusters obtained from the multivariate analysis. The Gamma function was used for the first level; this tool is very useful for ecological data with values above zero and kurtosis skewed to the right (^{Bolker 2008}), as is the case of population density. Differences in deer density for the different levels of analysis (subparagraphs 2, 3 and 4 above) were compared using a one-way analysis of variance and tested a posteriori with Tukey’s HSD multiple comparison test.

To determine the relationship between deer density and vegetation variables (Arb_den, Arb_alt, Arb_ab, N0, N1, H, E1, S, cob_50, and cob_100), we first examined whether there was a correlation between some of these variables using the R pairs function. This preliminary analysis allowed the elimination of six of these variables. Subsequently, we analyzed the relationship between variable-dependent density through a generalized linear model. All statistical and graphical analyses were conducted using the statistical program R 3.1.2 (^{R Development Core Team 2015}).

RESULTS

Classification of vegetation. A total of 1,756 trees of 101 species belonging to 49 botanical families were counted. Four groups or clusters of transects were obtained, which correspond to the locations studied (Figure 2a), as follows: group 1, Chiquero; group 2, Corralillo; group 3, Pensador; and group 4, Monte Grande; the Porvenir locality can be classified into groups 1 and 2 (Figure 2b). The similarity between transects of the same group and the difference between groups can be explained by the difference in the abundance of tree species. Group 1 was represented mainly by Acacia pennatula, Pithecellobium lanceolatum and Aldama dentatathe; group 2 was the most heterogeneous one, represented mainly by Malvaviscus arboreus, Tabebuia rosea, Croton glabellus and Croton sp; group 3 is represented mainly by Cedrela odorata. According to the degree of similarity, group 4 showed the lowest similarity relative to all other groups, and is represented by Guazuma ulmifolia, Croton glabellus and Ehretia tinifolia (Figure 2c).

Figure 2 Classification of transects using the Ward’s clustering method (a). Spatial distribution of the 40 transects classified into four groups (b). The colors in both figures match. Spatial distribution of the three most abundant tree species in the transects in each of the five study locations (C). Circle size is proportional to tree abundance. 

With regard to the comparison of vegetation types, those corresponding to agriculture and secondary medium forest vegetation were different vs. agriculture-pastureland and secondary vegetation with high-forest elements. The significant variables were tree height (F = 3.37, P = 0.029), basal area (F = 4.88, P = 0.006), tree richness (F = 5.01, P = 0.005), understory richness (F = 4.88, P = 0.006) and diversity of the tree stratum (F = 4.56, P = 0.008). On the other hand, at the level of groups and localities, group 3 corresponding to Pensador had the lowest density (F = 5.91, P = 0.002), height (F = 6.86, P = 0.001), basal area (F = 22.58, P = 0.002) and cover at 100 cm (F = 6.05, P = 0.001). Group 1 (Chiquero and Porvenir), group 2 (Corralillo and Porvenir) and group 4 (Monte Grande) did not differ with each other.

Relationship between deer density and vegetation. There were a total of 421 pellet groups in 39 of the 40 transects sampled in the five study localities. Average density for the entire region was estimated at 4.0 ± 2.8 ind/km². Density varied between localities, being statistically lower in Pensador (Figure 3a; F = 5.1, d .f. = 4, 35, P = 0.002). With regard to vegetation types, a trend was observed toward a lower deer density in agricultural areas and a higher density in forest patches (Figure 3b; F = 3.39, d. f. = 3, 36, P = 0.03). With regard to groups based on similarity in composition, the statistically lowest density was estimated for group 3 (Figure 3c; F = 4.85, d. f. = 3, 36, P = 0,006).

Figure 3 Estimates of the density (ind/km²) of the white-tailed deer in the study localities (a), vegetation types (b), and the groups (clusters) obtained in the multivariate analysis (c). Black lines inside boxes represent the value of the median. Lowercase letters above each bar represent the a-posteriori comparisons of the Tukey HSD test; same letters indicate non-different estimates. 

Of the ten vegetation variables originally calculated, six were eliminated for displaying a high colinearity. The coefficients of the first GLM suggest that tree density and Hill’s diversity index do not make a significant contribution (Table 3). The second model indicates that a higher percent understory cover and a lower basal area of the tree stratum predict a higher white-tailed deer density in the study area (Figure 4).

Figure 4 Relationship between understorey cover and basal area with white-tailed deer density. The red line represents the curve of best fit according to the linear models. The correlation coefficient (r) is shown. 

DISCUSSION

The white-tailed deer was recorded in 98 % of the transects sampled, with an estimated mean density of 4 ind/km² for the whole area. This density lies within the range reported for other regions of the state of Veracruz, i. e. 1 to 8 ind/km² (^{Contreras-Verteramo et al. 2014}; ^{Serna-Lagunes et al. 2013b}; ^{Mandujano and López-Tello 2015}). Deer density in the landscape studied is relatively similar to the one reported for other regions with secondary vegetation alternating with remnants of medium forest in southeastern Mexico (^{González-Marín et al. 2008}; ^{Bello et al. 2004}).

The use of generalized linear models enabled to predict a higher density of white-tailed deer in areas where the protection cover is high and the tree stratum is dominated by species with relatively low basal areas, i. e. thin trees such as Acacia spp., Mimosa spp., Cedrela odorata and Guazuma ulmifolia. Several of these genera and species are consumed by deer elsewhere (^{Villarreal-Espino 2006}; ^{Vázquez et al. 2015}). The five study localities were classified according to the composition and abundance of the main tree species. In particular, the highest estimated densities corresponded to Monte Grande, Porvenir, Chiquero and Corralillo, with a more diverse and dense vegetation composition and structure. In the landscape studied, the highest cover was found in patches with secondary vegetation resulting from slash-and-burn agriculture, which allows the regeneration of the tree vegetation after two or three years of use for agriculture (^{Reyes 2009}). In contrast, the lowest density was estimated for Pensador, a locality characterized by less favorable conditions for the white-tailed deer due to increased livestock raising and greater poaching pressure according to some persons. Other studies in forests have also documented that deer densities are higher where the vegetation provides this and other deer species with food and a dense cover to meet their needs for protection against harsh temperatures and predators (^{Ortiz-Martinez et al. 2005}; ^{Medina-Torres et al. 2008}; ^{Flores-Armillas et al. 2013}).

The white-tailed deer is exploited by local inhabitants, but as of today the number of animals hunted in the region has not been quantified, nor whether the rate exceeds a potential sustainable use. According to the experience in central Veracruz, the white-tailed deer could be managed under the diversified livestock scheme to contribute to conservation and use, incorporating other productive activities such as ecotourism, hunting tourism and scientific tourism, as an alternative to generate income for UMA owners, contributing to the improvement and conservation of ecosystems at the regional level (^{Serna-Lagunes et al. 2013b}). According to ^{Mandujano (2010)}, the sustainable development model in extensive UMAs in tropical regions is feasible only if these incorporate areas of several thousands of hectares at the landscape level (^{Mandujano and González-Zamora 2009}). In this sense, a recent experience in central Veracruz suggests that the management of this deer species is not sustainable in UMAs of less than 3000 hectares (^{Mandujano and López-Tello 2015}). Another relevant aspect is that intense and constant work with human communities is required in order to incorporate them to the management of this animal resource (^{Lopez-Tellez et al. 2016}). In addition, considering that the study region is still home to human communities of different ethnic origin, conducting intense work from an ethno-zoological perspective is suggested (^{Mandujano and Rico-Gray 1991}; ^{Naranjo et al. 2004}).

Similar to what has happened in other regions of the state of Veracruz (^{Landeros-Sánchez et al. 2011}), the landscape studied is currently made up of a complex mosaic of temporary and permanent agriculture, cultivated pastureland for bovine cattle raising, and remnants of secondary vegetation with the tree stratum of medium semi-deciduous and high tropical forest. In this sense, and considering that sampling was performed only during the dry season of 2014, a very interesting aspect to address in future studies refers to the structural complexity of landscape and the presence of relatively large patches of tree vegetation, which would benefit the white-tailed deer and other wildlife species, similar to what has been documented for other regions of southeastern Mexico (^{García-Marmolejo et al. 2015}; ^{Contreras-Moreno et al. 2016}).