Data analysis
It was difficult to calculate the burrow density, as the area is a mosaic of different geological forms (Kabala et al., 2021). Some parts of the valley could be unsuitable for burrow forming, or avoided by marmots for other reasons, affecting burrow density. Therefore, we had to determine the area for burrow density measurement arbitrarily. We also calculated the mean distance from each burrow to the nearest burrow and from the five nearest burrows.
We excluded burrows that were located on the edge of the highest river terrace from spatial analysis of aerial images. For the remaining 60 burrows, we analysed land cover in round buffers with a radius of 20 m around the burrow entrance. We divided the buffer area of the selected burrows into hexagons with a short diagonal of 25 cm, corresponding to an area of 0.054 m2. We assigned each hexagon to one of four categories: burrow entrance, burrow mound, bare soil or vegetation. The categories reflected the dominant type of cover (>50% area) within the hexagon. There was only one burrow entrance hexagon per burrow. Burrow mounds were easily distinguishable due to their grey colour and determined through visual inspection of the orthophotomosaic. Vegetation was identified using raster analysis with the Red-Green-Blue Vegetation Index (Bendig et al., 2015). All hexagons that did not fall in the three mentioned categories were described as bare soil. Distance measurements were done between centroids of hexagons. For each buffer around a burrow, we measured the following parameters: (1) total vegetation cover of the area free of mounds as: vegetation hexagons/(vegetation hexagons + bare soil hexagons) (2) mound area (3) mean and standard deviation of the distance of soil and vegetation hexagons from the burrow entrance. For the analysis of the relationship between vegetation cover and the distance from the nearest burrow entrance, we selected burrows with a total vegetation cover greater than 2%. From the remaining 47 burrows, each was analysed separately. In each buffer around a burrow, we divided hexagons into 20 equal intervals based on their distance from the nearest burrow entrance, each interval 1 m wide. Due to the small number of hexagons in the first interval, we merged it with the second interval, giving a total of 19 intervals. In each interval, we calculated the vegetation cover of the area free of mounds. We identified outliers based on the z-score, which was calculated as: z = (value – mean value from the buffer)/SD from the buffer. Data points with a z-score below -3 or above 3 were excluded from further analysis. We used a binomial linear model to calculate the adjusted R2 value for the relationship and plot the data. Having four values for every buffer (total vegetation cover, the mean distance of vegetation hexagons from the nearest burrow entrance, mound area and R2adjusted for the relationship between the vegetation cover and the distance from the nearest burrow entrance), we calculated the Pearson correlation coefficient (r) between each pair of parameters.
Before statistical analysis of the results of chemical analyses, we checked whether the data met the requirements of parametric tests: normal distribution (using the Shapiro-Wilk test) and equality of variances (using the F-test for two groups and the Bartlett test for multiple groups). We performed the Kruskal-Wallis test to compare the chemical parameters of plant species and the post-hoc Dunn’s all-pairs test with Holm adjustment of p-values. We used the two-sidedt -test (for normally distributed data with equal variances), the Welch’s t -test (for normally distributed data with unequal variances) or the Wilcoxon signed-rank test (for non-normally distributed data) to check for differences in the abovementioned parameters between plants of the same species from the high and low cluster, and between faeces from the high and low cluster. We calculated the Spearman’s rank correlation coefficient (rho) between the measured parameters in biomass and the distance of the plant from the burrow entrance, separately for each plant species. Based on N and P content in plant biomass, we calculated the N:P ratio.
We used QGIS v. 3.10.14 software for spatial analysis (QGIS Development Team) and map preparation, and R v.4.0.5 software (R Core Team 2021) for statistical analysis, using the build-in ‘stats’ package. We used the ‘ggplot2’ package (Wickham, 2016) to prepare graphs.