Effects of sex and season (breeding and non-breeding) on microhabitat selection in Stejneger’s bamboo pitviper (Viridovipera stejnegeri)

Study area

Data were collected in Huangshan City, Anhui Province, eastern China (117°23′–118°55′E, 29°24′–30°24′N; Fig. 1), the study area covers 25 square kilometers. The region is characterized by a subtropical monsoon climate with abundant heat and moisture. Annual average temperature ranges from 15.5 to 16.4 °C, and average annual precipitation ranges from 1,395 to 1,702 mm, predominantly falling from May to August (summer), with reduced rainfall from September to November (autumn), although it does not reach an extreme level of aridity (The People’s Government of Anhui Province, 2023). The vegetation is primarily composed of subtropical evergreen broad-leaved forests, the understory is primarily composed of perennial herbs (Poaceae and Urticaceae) and ferns, with minimal variation in vegetation diversity between seasons.

Data acquisition and variable assignment of categories

Data were collected during two periods, June–July 2023 and September–October 2023, from approximately 20:00 to 24:00. We extensively searched for V. stejnegeri in the survey area, and manually collected snakes. Each located snake was injected with an independent electronic tag (EM4305; wnLikes, Shenzhen, China), and its unique identification number, latitude and longitude coordinates, and sex were recorded. Previous studies on the closely related species Trimeresurus macrops, have revealed that this group typically exhibits a relatively small hunting range and home range. Therefore, reference to T. macrops, at each sampling location, a 4 × 4 m quadrat was established (Strine et al., 2018). Within these quadrats, 13 habitat factors were measured, with detailed descriptions and measurement methods provided in Appendix S1. A total of 92 snakes were collected (females = 34, males = 58), all snakes were released after data collection. The data were obtained only at the time of the first collection.

Sex identification was performed via cloacal probing by gently pressing the region located five cm below the vent at the tail towards the anterior direction of the cloaca using the thumb. The detection of an everted hemipenis was used to determine the sex of the individual. Existing research indicated that V. stejnegeri typically given birth in August (Huang, 1980; Huang et al., 2015). Therefore, data collected between June and July were designated as breeding season data, while data from September to October were considered non-breeding season data. During each survey period, ArcGIS v10.8 was used to randomly generate 50 pseudo-absence points within the survey area, excluding unsuitable habitats for V. stejnegeri such as highways, open water bodies, and cliffs, as well as points with the presence of V. stejnegeri within a 50-meter radius (Mizsei et al., 2024). The remaining points served as pseudo-absence points for measuring the aforementioned habitat factors (Zhang, 2020). All animal procedures were carried out in accordance with and approved by the Animal Care and Use Committee at Yibin University (animal protocol number: YBU2020007). We also adhered to the ARROW (Animals in Research: Reporting On Wildlife) guidelines.

Data processing and analysis

The collected data were categorized into four groups based on season and sex: breeding season males (BM), breeding season females (BF), non-breeding season males (NBM), and non-breeding season females (NBF).

Before starting data analysis, we checked and removed all extreme outliers and used the chain equation multivariate interpolation method to impute missing values. And all data were logarithmically processed before analysis to improve normality. All data were tested for normality and homogeneity of variance.

To analyze the habitat selection preferences of V. stejnegeri for each habitat factor, the Vanderploeg (W_i) and Scavia (E_i) selection indices were calculated (Vanderploeg & Scavia, 1979):

${\displaystyle W_i=\left(\frac{r_i}{p_i}\right)/\sum \left(\frac{r_i}{p_i}\right)}$ ${\displaystyle E_i=\left(W_i-\frac{1}{n}\right)/\left(W_i+\frac{1}{n}\right)}$

where r_i denotes the number of quadrats selected by V. stejnegeri with characteristics i; p_i denotes the total number of quadrats in the environment with characteristics i; and n refers to the level or category number of a particular habitat factor. E_i ranges from −1 to 1, with E_i = 1 indicating strong preference, 0.1 < E_i < 1 indicating selection, −0.1 < E_i < 0.1 indicating random selection, −1 < E_i < −0.1 indicating avoidance, and E_i = −1 indicating no selection.

To analyze the importance of habitat factors for each group, factor autocorrelation analysis was performed to exclude highly correlated factors (see Appendix S3). Random forest models were then constructed using the habitat data for each group as well as the pseudo-absence points data. We used viper presence/pseudo-absence as the dependent variable and employed the R package “randomForest” to construct random forest models, employed the mean decrease Gini index to evaluate the importance of each factor (Zhang et al., 2020).

Multiple linear models were used to analyze the interactive effects of season and sex on various factors. For factors showing significant effects, post hoc multiple comparisons (Least significant difference, LSD) were conducted among the four groups, while for factors without significant effects, differential analyses were restricted to within either seasonal or sexual groups. Prior to within-group analyses, tests for normality and homogeneity of variance were conducted on continuous variables. Continuous variables that met the assumptions of normality and homogeneity of variance were analyzed using one-way analysis of variance (ANOVA), while those that did not meet these criteria were analyzed using the Mann–Whitney U test. Discrete variables were examined using chi-square tests to identify differences within groups. We also calculated the mean and standard deviation within each group. Furthermore, partial dependence plots were employed to visually depict the relationship between habitat factors and predicted probability of occurrence for each group, as determined from the random forest models (Hastie, Tibshirani & Friedman, 2005).

To analyze overall differences and associations in microhabitat selection among groups, the microhabitat niche width Levins index (B) for each group was calculated using the spaa package in R (R Development Core Team, 2018): ${\displaystyle B=\frac{1}{\sum _{j=1}^R{\left(P_j\right)}^2}}$

where P_j represents the proportion of individuals found at sampling point j and R denotes the total number of sampling points (Simpson, 1949). Subsequently, microhabitat niche overlap between groups was quantified using the Levins index (O_ik): ${\displaystyle O_{ik}=\frac{\sum _{j=1}^R{P}_{ij}{P}_{kj}}{\sum _{j=1}^R{\left(P_{ij}\right)}^2}}$

where P_ij and P_kj represent the abundance of groups i and k at sampling point j, respectively, and R is the total number of sampling points (Levins, 1968).

To calculate the significance of niche overlap, we randomly resampled data regardless of group identity sort them to two group, calculate niche overlap and compare observed value by a t-test.

Principal component analysis (PCA) was conducted to determine the eigenvalues of each principal component (PC). PCA scatter plots were created using the PCs with the highest eigenvalues as axes. The distribution of points within 95% confidence ellipses for different groups on the scatter plot was analyzed to assess whether significant differences existed among the groups. To further validate the results, permutational multivariate analysis of variance (PERMANOVA) (ADONIS) was conducted using the vegan package in R, with the degree of overlap among groups in the graphical representations assessed to determine whether significant differences existed (R Development Core Team, 2018).

All analyses were conducted with a confidence level of 0.05. Data analysis was performed using R v4.2.3, and graphs were plotted using R v4.2.3 and Origin 2022 (R Development Core Team, 2018).

Microhabitat preferences and factor importance among groups

A total of 62 presence points were collected during the breeding season (males = 43, females = 19) and 30 during the non-breeding season (males = 15, females = 15), all snakes are adults (SVL: 57.04 ± 5.19 cm).

The Vanderploeg (W_i) and Scavia (E_i) selection indices revealed no significant differences among the four groups in their preferences for factors such as temperature, humidity, vegetation height, vegetation coverage, slope, distance from water, distance from residential sites, landscape habitat, and vegetation type. However, compared to other the groups, NBF displayed a preference for lower altitudes and slopes, and for locations furthest from roads. Additionally, all groups showed either no preference or extremely weak preference for “aspect” factor, leading to its exclusion in subsequent analyses (Fig. 2).

The mean decrease Gini index identified the key factors influencing habitat selection. For the BM group, important factors included humidity, temperature, distance from water, and distance from roads. The BF group was influenced by humidity, temperature, and distance from water, while the NBM and NBF groups were influenced by humidity and distance from water (Fig. 3).

Seasonal (breeding and non-breeding) and sexual differences

Multiple linear model analyses revealed that, apart from altitude, there were no significant interactions between sex and season for other factors (Appendix S2).

Mann–Whitney U and chi-square tests demonstrated significant differences between males and females only in distance from residential sites (M vs F: 414.98 ± 170.96 vs 338.65 ± 148.48; P = 0.026) and slope position (F vs M: mid-slope position vs down-slope position; P = 0.004; Table 1). Seasonal differences analyses only indicated a significant difference in temperature preference (M vs F: 21.95 ± 0.81 vs 23.71 ± 0.84; P < 0.001; Table 1), with no significant differences in preferences for other factors. Specific differences among groups are illustrated in Fig. 4.

Multiple comparisons analysis revealed that the altitudes selected by NBF differed significantly from those selected by NBM (P < 0.001), BM (P = 0.001), and BF (P = 0.009). However, there were no significant differences in altitude selection among NBM, BM, and BF. Partial dependence plots indicated that NBF selected lower altitudes than the other groups, with a narrower range of altitude preferences (NBF: 190.73 ± 19.06, BM: 212.09 ± 15.91, NBM: 211.42 ± 17.63, BF: 208.40 ± 22.77, Fig. 4A).

Overall differences in microhabitat selection among different groups

Microhabitat niche analysis indicated that the microhabitat niche widths varied among the different groups, with NBM exhibiting the largest (2.935) and NBF the smallest (1.895). As shown in Table 2, the greatest microhabitat niche overlap was found between NBM and BF (98.63%), while the smallest overlap was found between NBF and BM (53.86%). The results of the paired sample t-test between randomly resampling and observed values show that the niche overlap index of the observed values is slightly lower than that of randomly resampling (79.74 ± 18.60 vs 85.61 ± 14.89), but the difference is not significant (t = 0.533, P = 0.617). This indicates that niche differentiation among the groups is not pronounced.

The PCA plot indicated no significant differences among the four groups (Fig. 5), supported by ADONIS analysis (Fig. 6).