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.