Introduction

Co-occurrence processes between different species depend on the spatial scale at which they perceive the habitat (^{Morris 1987}; ^{Barrio and Hik 2013}). For instance, it has been proposed that small mammals possibly perceive the habitat structure/heterogeneity at smaller scales than medium-sized or large mammals (i. e., microhabitat; ^{Morris 1984}, ¹⁹⁸⁷; ^{Chesson 2000}; ^{Whittaker et al. 2001}). This is relevant since spatial scale may alter species assemblage’s patterns perception and the order of importance of the explanatory variables of these patterns (^{Morris 1987}; ^{Whittaker et al. 2001}). Likewise, coexistence requires species to be different in the way they affect and are affected by competitors and available resources, resulting in niche differences or average fitness differences between species (^{Chesson 2000}; ^{Chen et al. 2020}). Indeed, differences in the ecological niches occupied by the species within a community act to stabilize the system, with larger differences promoting coexistence ^{(Chen et al. 2020}). Therefore, determination of the mechanisms for coexistence among several species within a community is of basic ecological interest.

Habitat heterogeneity (i. e., different biotic and abiotic components that shape the habitat architecture) is expected to increase species coexistence of small mammals, as they increase the number of microhabitats that may be occupied by species with different environmental requirements (^{Rosenzweig and Winakur 1969}; ^{Cramer and Willig 2002}; ^{Corbalán and Ojeda 2004}; ^{Schreiber and Killingback 2013}; ^{Novillo et al. 2017}). Likewise, more available microhabitats offer more shelters for mice to hide from predators, and also provide more diversity of food resources (^{Corbalán and Ojeda 2004}; ^{Traba et al. 2010}; Novillo et al. 2017). However, temporal changes in the habitat structure have also a substantial impact on species coexistence (^{Valladares et al. 2015}). For example, changes in the availability of food resources and refuges for species through space and time, affect its population size and how they interact (^{Valladares et al. 2015}). Indeed, temporal variation in habitat heterogeneity should increase available niche spaces, allowing more species to coexist (^{Currie 1991}). It has been reported that temporal fluctuations in habitat structure can stabilize species coexistence via the “storage effect” (^{Chesson 2000}), when inter and intra-annual variation in climate or resource availability favors one group of species over others (^{Zavaleta et al. 2003}). Similarly, it has reported that the sex of individuals could contribute to the coexistence of two competitive species (^{Chesson 2000}). For instance, it has found a positive effect on coexistence when differences in competitive ability among conspecific individuals (e. g., competition for resources between males vs. males or males vs. females of the same species) can break down competitive hierarchies of species, such that intraspecific competition is stronger than interspecific competition (^{Chesson 2000}; ^{Hubbell 2005}; ^{Fridley et al. 2007}; ^{Uriarte and Menge 2018}).

In this study, we focus on whether different microhabitat elements and their spatial and temporal changes shape the frequency of two syntopic (animals that may use the same habitat; ^{Hart et al. 2018}) and congeneric species of deer mice (P. difficilis and P. melanotis). These two model systems were selected for several reasons. First, these two congenerics co-occur in similar habitats and environmental conditions along its distribution (^{Álvarez-Castañeda 2005}; ^{Fernández et al. 2010}). In the Transmexican Neovolcanic Belt, they cohabit in the understory of mixed and coniferous temperate forests, such as in the Desierto de los Leones National Park (DLNP; ^{Castro-Campillo et al. 2008}), where this study was conducted. Second, since the DLNP is located at the edge of an ever-growing megalopolis (Mexico City), it is important to understand what elements of the habitat shape the frequency and coexistence of these kind of species to elaborate better strategies to reduce the human impact on natural areas and their wild inhabitants. Urban growth produces fragmentation of natural microhabitats that are important to small mammals, which in turn play a fundamental role in the dispersal of seeds and as habitat architects within the forest.

Third, these two congeneric species have different body sizes (P. diffilis is larger than P. melanotis (see^{Álvarez-Castañeda 2005}; ^{Fernández et al. 2010}). This morphological differentiation represents an opportunity to compare their different requirements in relation to their respective ecological niche. Four, the reproductive season of P. difficilis occurs mainly during the dry season, while that of P. melanotis occurs in the rainy season (^{Castro-Campillo et al. 2012}; ^{De-la-Cruz et al. 2019}; ^{Salame-Méndez et al. 2018}, ²⁰¹⁹, ²⁰²⁰). Thus, it is possible that spatial and temporal variation in habitat heterogeneity, such as availability of food resources and shelters, could be related to the different reproductive seasons, relaxing its interspecific competition and promoting its coexistence (^{Kaufman and Kaufman 1989}; ^{Chesson 2000}). Nevertheless, while several studies have quantified variation in resource use (e. g., food, shelters, water) in heterogeneous environments among mice from the same genus (^{Kaufman and Kaufman 1989}; ^{Kalcounis-Rüppell and Millar 2002}; ^{Villanueva-Hernández et al. 2017}), there is still a lack of information about how temporal changes in the habitat structure/heterogeneity affect the frequency of congeneric mice species inhabiting at the same place (but see^{Kalcounis-Rüppell and Millar 2002}; ^{Hart et al. 2018}). Here, we measured different variables that compose the structure of the microhabitat and their changes during two rainy seasons (dry and rains), and we related this information with total and sex frequencies of two species of small mammals, P. difficilis and P. melanotis, that co-occurs in a temperate forest of Central Mexico.

Materials and methods

Study area. The study area (19º18’17”N, 99º19’14” W at 2,289 masl) is located in a mixed temperate forest of coniferous and broad-leaved trees at Desierto de Los Leones National Park (DLNP) in Mexico City (^{CONANP 2006}). This forest is part of the Trans-Mexican Neovolcanic Belt (^{CONANP 2006}). The rainy season occurs from summer through early fall (June to October) with a monthly average precipitation of 252.92 ± 28.01 mm, and an average monthly temperature of 11.72 ± 0.53 °C (^{CONANP 2006}). In contrast, the dry season occurs from fall through winter (October to February); the monthly average precipitation is 13.2 ± 3.11 mm, and the average monthly temperature is 8.97 ± 0.68 (^{CONANP 2006}).

Habitat heterogeneity during dry and rainy seasons. To assess how temporal changes in the habitat heterogeneity affect the frequency of capture of each Peromyscus species, we set a plot of 9,900 m² (110 x 90 m; Figures 1a, b). Within this plot, we placed 120 sampling stations every ten meters, along 12 columns (A to L) and ten rows (1 to 10; Figure 1a). Each sampling station was marked with buried wooden stakes. In each sampling station, we delimited an “influence zone” of 2.5 m² (Figure 1b). Within these influence zones of every sampling station, eight fixed transects were set in a clockwise conformation to measure seven variables that qualify as components of the vertical and horizontal structure of the habitat (^{Morris 1984}; ^{Jorgensen 2004}; ^{Villanueva-Hernández et al. 2017}). These variables are indicators of possible shelters from predators, spaces for resting and mating, and food resources (^{Jorgensen 2004}). We applied the Canfield’s Line Intercept (CLI) method (^{Canfield 1941}) in each transect of the influence zone to measure the percentage of vegetation coverage at three different heights (10, 35, and 100 cm; VC10, VC35, and VC100, respectively). We also counted the number of all herbaceous plants (H) and the number of all woody plants (W). Likewise, plant species richness (SR) was registered as the number of different plant species found within each influence zone (Figure 1b). The percentage of logs (Logs) covering the ground surface (fallen trees) more than one meter long and ten centimeters in diameter was also registered using Canfield’s method. The CLI method is based on the measurement of all plants and objects intercepted by a transect, and the length/coverage of each plant or object that is touched by the line is registered. All habitat features were sampled once during the most representative months of the rainy (July 2017) and dry seasons (February 2018; ^{CONANP 2006}).

Figure 1 a) Location of the sampling plot in a temperate mixed forest at the Desierto de los Leones National Park, Mexico City. b) Configuration of the experimental plot area (grid) where 120 sampling stations were set with its corresponding influence zone of 2.5 m². Within each influence zone, eight fixed transects were set in a clockwise conformation (North to Southwest) to measure the variables related to the habitat heterogeneity. Green triangle = center of the sampling station where each Sherman trap was also set. 

Mice trapping. The Peromyscus mice were captured alive for ten months to include data for the dry (October 2017 to February 2018) and rainy (March 2018 to July 2018) seasons. We set a single live trap (H. B. Sherman, Inc., Tallahassee, FL 32303, USA), baited with oat flakes and vanilla extract at each sampling station of the plot (n = 120). Traps were set for two consecutive nights each month with a total of 20 capture events (10 months x 2 nights; 2,400 night/traps). Mice were marked in their abdomen with gentian violet to avoid overestimation of capture frequencies. Despite this temporal mark did not allow us to identify the mice at individual level (e. g., as with an id code to trace the movements or range home of a certain individual), it turned out to be a practical way to calculate mice abundance accurately by pulling out the recaptured individuals. At the end of the samplings, we could have an estimate of the mice abundance for the dry and rainy seasons. In addition, to avoid recurrence behavior (e. g., mice returning to the traps for bait) or shyness (e. g., mice avoiding traps due to other mice odors), all the traps were thoroughly cleaned and randomly oriented within each sampling station in each capture event. Handling of the mice was following the guidelines of the American Society of Mammalogists (^{Sikes 2016}). To prevent hypothermia during capture, we placed 3 to 5 cotton balls inside the trap and put the traps inside open plastic bags. We recorded the species and sex of every mouse (^{De-la-Cruz et al. 2019}). If a mouse died overnight while being trapped (n = 3), it was skinned and prepared as a study specimen, and incorporated as voucher specimen in the Mammal Collection of the Universidad Autónoma Metropolitana-Iztapalapa. A scientific collecting permit, SEMARNAT-08-049-B, was issued to Alondra Castro-Campillo (ACC) by DGVS, SGPA-09712/13, SEMARNAT, MEXICO.

Statistical analysis. All statistical analyses were performed using the JMP statistical package (v. 14.0; SAS Institute). Plotting was made using ggplot2 (^{Wickham 2016}) in RStudio version 1.1.463 (^{R Core Team 2020}).

Differentiation of habitat heterogeneity between dry and rainy seasons. Prior to analyses, all variables (habitat elements and mice captures) were log-transformed (log n + 1) to meet normality assumptions. A Student’s t-test was used to evaluate mean differences between the rainy and dry seasons for each of the seven variables of the habitat heterogeneity (see above).

Relationship between habitat heterogeneity and frequency of each Peromyscus species during the dry and rainy seasons. Prior to analyses, variables representing the habitat’s heterogeneity and structure were standardized to a mean of zero and a standard deviation of one (X^ = 0, SD = 1). Generalized linear models (GLMs) were used to evaluate the relationships between habitat heterogeneity and the capture frequency of each Peromyscus species. All GLMs were performed with 1,000 iterations of the maximum likelihood method. The GLMs described hereafter were selected based on the statistical significance of the model and the lowest corrected AIC values, that is, models that best explained the relationship between the variables (^{Akaike 1974}).

To evaluate the relationship between the habitat heterogeneity and the frequency of each Peromyscus, two GLMs (link = logarithmic, distribution = Poisson) were constructed using the respective abundances of P. difficilis and P. melanotis as response variables, whereas the vegetation coverage at 10 and 100 cm (VC10 and VC100, respectively), the number of all herbaceous plants (H), the number of all woody plants (W), plant species richness (SR), the percentage of logs (Logs), season, and their interaction were used as predictors. Adding the interaction between season and the covariates in the models, allowed us to assess whether the effect of the habitat heterogeneity or structure differed depending on the season of testing (cf. ^{Zar 1999}). A similar GLM was carried out using the abundance of females or males of each Peromyscus as response variables.

Since habitat heterogeneity depends on the interaction between different biotic and abiotic elements of the habitat, we also carried out principal component analyses (PCA) for each season, using the seven habitat indicators (variables) to produce new functions that could explain the microhabitat heterogeneity in a more detailed way. Likewise, the PCAs also helped to reduce data dimensionality. One PCA was performed for each season. Since in the dry season, the first two components explained 61.23 % of the variance and in the rainy season, the two first components explained 73.88 % of the variance, we used only these two principal components for subsequent analyses (see below).

A GLM (link = logarithmic, distribution = Poisson) was constructed where the response variable was the abundance of P. difficilis and/or P. melanotis and the two principal components, season and their interaction as covariate effects. A similar GLM was carried out using the abundance of females or males in each Peromyscus species as response variables. The generalized linear coefficients (viz. βi; ^{Lande and Arnold 1983}) obtained from the GLMs represent the strength and direction of the relationships acting directly on the frequency of the species in comparable units (standard deviations).

Results

Frequency of captures. We captured a total of 516 mice of both Peromyscus species for all the study. During the dry season, the total number of captured individuals was 312: 174 individuals for P. difficilis (120 males and 54 females) and 138 P. melanotis individuals (66, 72). During the rainy season, the total number of captures was 204; 114 individuals of P. difficilis (80, 34), and 92 individuals of P. melanotis (58, 34).

Changes in habitat structure between the dry and rainy season. The mean difference between seasons was significant for most of the habitat features (Table 1). The mean of almost all habitat variables was higher in the rainy season. Only the coverage of wood logs on the ground surface between seasons did not change (Table 1, Figure 2).

Relationships between frequency of capture for the Peromyscus species and the seven variables of habitat heterogeneity/structure. The GLM between the seven variables that describe the microhabitat heterogeneity/structure and the frequency of P. difficilis was significant (L-R chi-square₁₅ = 37.16, AICc = 451.30, P = 0.0012, Table 2, Appendix 1). However, only the percentage of logs on the ground (positive effect) and season (positive effect) as predictors were significant. The GLM between the seven habitat variables and P. difficilis females’ frequency was not significant (L-R chi-square₁₅ = 15.56, AICc = 263.26, P = 0.4118, Appendix 1, 2). In this model, the percentage of herbaceous plants (H) showed a significant negative effect on the frequency of females of P. difficilis. However, this result must be interpreted with caution. In contrast, the effect of the seven habitat variables on the frequency of males of P. difficilis was significant (L-R chi-square₁₅ = 44, AICc = 365.51, P = 0.0001, Appendix 1, 2), being percentage of logs covering the ground and the interaction season × woody plants significant predictors (positive effects).

Figure 2 Median differentiation of the seven variables that describe the habitat heterogeneity between dry and rainy seasons. Dots represent each sampling station. Y-axis is in log scale. 

The GLM between the seven variables and the frequency of P. melanotis was significant (L-R chi-square₁₅ = 36.37, AICc = 401.37, P = 0.0016, Table 2, Appendix 1). Significant effects included the season and plant species richness. Both effects were positively related to the frequency of P. melanotis. The GLM between the seven habitat variables and the frequency of P. melanotis females was nearly significant (L-R chi-square₁₅ = 24.93, AICc = 277.54, P = 0.0508, Appendix 1, 2). In this model, only the season was significant. The effect of the seven habitat variables on the frequency of P. melanotis males was significant (L-R chi-square₁₅ = 34.59, AICc = 295.11, P = 0.0028, Appendix 1, 2). In this GLM, the effect of the species richness and the interaction between season and species richness were significant and positively related to P. melanotis males. On the other hand, the number of woody plants was significant and negatively related to P. melanotis males.

Relationships between capture frequency and principal components. Principal component analysis for the dry season showed that the first two components explained almost all the variance (61.23 %). The highest loadings in the first component (PC1) were variables related to vegetation including plant species richness. Hence, this new function describes those sampling stations with high vegetation coverage and species richness. In contrast, the principal component two (PC2) was only related to the coverage of logs on the ground, representing the fixed and stable elements of the habitat. In the rainy season, the first two principal components explained 73.88 % of the variance. As in the dry season, PC1 was related to vegetation and species richness and PC2 with logs’ coverage (Appendix 3).

The GLM between the frequency of P. difficilis and the principal components was significant (L-R chi-square₅ = 30.01, AICc = 436.37; P = 0.0001; Table 3, Figure 3 a, b). However, significant predictors only included the PC2, indicating a positive relationship between the frequency of P. difficilis and logs’ coverage during the two rainy seasons. The GLM between the frequency of P. difficilis females and the principal components was not significant (L-R chi-square₅ = 3.15, AICc = 253.60; P = 0.6767) (Appendix 2). The effect of the principal components on the frequency of P. difficilis males was significant (L-R chi-square₅ = 36.28, AICc = 351.15; P = 0.0001). However, the only significant and positive effect was the PC2 (logs on the ground; Appendix 2).

The GLM between the abundance of P. melanotis and the principal components was significant (L-R chi-square₅ = 22.46, AICc = 393.20; P = 0.0004) (Table 3, Figure 3 c, d). The significant and positive predictors were the PC1 (vegetation) and season in this model. The GLM between the frequency of P. melanotis females and the principal components was significant (L-R chi-square₅ =12.93, AICc = 267.46; P = 0.0240; Appendix 2). The significant effects also included the PC1 and season (positive effects). The GLM between the frequency of P. melanotis males and the principal components was also significant (L-R chi-square₅ =17.92, AICc = 289.70; P = 0.0030; Appendix 2). The significant effects included the PC1 (positive effect) and PC2 (negative effect).

Figure 3 Relationships between the abundance of Peromyscus difficilis or Peromyscus melanotis with the principal components. p-values are showed in the plots. 

Discussion

Our results revealed different relationships between the habitat heterogeneity and the frequency of both Peromyscus species. Likewise, the habitat heterogeneity/structure had a different effect on male and female frequencies of both species of Peromyscus. First, while P. difficilis was more frequently captured in sampling stations with high presence and coverage of logs on the ground, individuals of P. melanotis were usually captured in sampling stations with high vegetation coverage and plant species richness. It has been found that spatial heterogeneity in habitat quality reverses the competition between two sympatric species (^{Chen et al. 2020}). Thus, it seems that both Peromyscus species have different habitat requirements and preferences that could be facilitating their coexistence in the study zone (^{Pianka 1973}; ^{Chen et al. 2020}).

Why P. difficilis had a higher association in places with a high presence of logs? Logs on the ground surface represent small patches of microhabitat with food sources, burrows and refuges for mice (^{Bellows et al. 2001}; ^{Bowman et al. 2015}; ^{Grelle 2003}; ^{Dalmagro and Vieira 2005}). Individuals of P. difficilis could also use the large fallen logs as safe pathways for fast and straight locomotion within the forest (^{Bellows et al. 2001}; ^{Grelle 2003}; ^{Dalmagro and Vieira 2005}). Indeed, fallen logs promote the structural heterogeneity of forests and may enhance positive interactions among species of small mammals (^{Bowman et al. 2015}). The relationship between P. difficilis and microhabitats with a higher number and cover of logs could also be related with the morphology and climbing habits of this species (^{Fernández et al. 2010}). For instance, the long tail of P. difficilis enables it to rush and climb along shrubs, trees, or logs (^{Bowman et al. 2015}), hence, increasing its preference for habitats with fixed elements (e. g., fallen logs), where mice can escape from predators or use holes in logs as burrows (^{Bowman et al. 2015}). In contrast, it is possible that P. melanotis - the species with smaller body size and more cursorial locomotion - prefers zones with higher vegetation coverage and plant species richness as strategy to avoid predation by aerial hunters or other predators, as well as to obtain food resources more quickly, since seeds may be concentrated under shrub canopies (^{Bowman et al. 2015}; ^{Grelle 2003}; ^{Dalmagro and Vieira 2005}). Likewise, we have found evidence that P. melanotis has a higher intermediate metabolism rate than P. difficilis (Salame-Méndez and Castro-Campillo, unpublished results). Thus, it is possible that due to its higher metabolic requirements, P. melanotis needs to disperse more frequently in search for food resources richer in sugar in areas with higher diversity and richness of plants (such as we observed in this study) that supply its higher metabolic requirements and promoting its ecological niche differentiation (^{Suarez and Welch 2017}; ^{Chen et al. 2020}).

Our results also revealed that temporal changes in the habitat heterogeneity only seem to affect the frequency of males of both Peromyscus species in the study zone, the sex more frequently related to territorial endeavors in mammals (^{Ostfeld 1990}). For instance, males of P. difficilis were positively associated with woody plants in the dry and rainy seasons. However, this association was stronger during the dry season. In contrast, males of P. melanotis were negatively associated with woody plants in the dry season and positively associated with plant species richness in both seasons, but the association was also stronger during the rainy season. Thus, it seems that there is a different microhabitat use by males of both Peromyscus species. Likewise, it is possible that the stronger association between woody plants and males of P. difficilis during the dry season could be related to the breeding season of this species. It has been reported that the breeding of P. difficilis occurs primarily during the dry season (^{De-la-Cruz et al. 2019}; ^{Salame-Mendez et al. 2020}). Hence, since males could be searching for females to breed during the dry season, it is possible that they are more easily captured than the females or that the density of males in the area increased due to breeding season. This same pattern could also explain the stronger association between males of P. melanotis and plant species richness or areas with higher vegetation coverage during the rainy season, where the breeding of P. melanotis occurs (^{De-la-Cruz et al. 2019}; ^{Salame-Mendez et al. 2020}).

Our findings also revealed that females of both Peromyscus species did not show a significant association with the microhabitat variables. One explanation is that females may be associated with other habitat variables not measured in this study. Likewise, it is possible that females spend more time in their burrows than males, and they only go out to obtain food and water for short periods, being more cautious and not necessarily associated with a specific habitat component. Indeed, it has been reported for several small mammal species that females have shorter home ranges than males (e. g. ^{Ribble et al. 2002}; ^{Flores-Manzanero et al. 2019}). Nevertheless, the lack of a significant association between females and habitat variables could also be related to a statistical bias (^{Zar 1999}), since we captured more males than females during both seasons, except for P. melanotis in the dry season where we captured six more females than males of this species.

Finally, this study provides evidence of how two congeneric mice species are affected by the habitat heterogeneity. Interestingly, although these two coexisting Peromyscus species can move along different microhabitats (i. e., logs, vegetation coverage), the fact that each one was associated with a different habitat component, suggests a different microhabitat use and selection (^{Chesson 2000}). Habitat partitioning is considered an important mechanism for coexistence among small mammal species (^{Schoener 1974}; ^{Stevens and Tello 2009}). Therefore, partitioning of some of the available resources (e. g. space, food, refuges, and time) could be enhancing the coexistence of these two syntopic Peromyscus species with similar environmental requirements (^{Shenbrot 1992}; ^{Chesson 2000}; ^{Novillo et al. 2017}). Overall, we demonstrate the importance of the study of habitat heterogeneity on the interaction of two syntopic species, offering some new insights into the mechanisms of spatial and temporal coexistence of two mice species at local scale.