The research was conducted at the Bakuriani ski resort, located 183 km west of Tbilisi (41.75°N, 43.53°E) on the slopes of the Lesser Caucasus, Georgia. The ski runs of the Didveli piste complex extend from 1,884 to 2,685 m a.s.l. across montane, subalpine, and alpine vegetation zones (Fig. 1).

Alpine skiing has been practised on this mountain since 1932; prior to that, the terrain was used as summer pasture. The grazing season lasts three months in the summer, with shepherds typically transferring their livestock to the mountainous regions from late May until September (REC Caucasus 2019).

The principal environmental pressures on soil and vegetation at present are grazing (cattle and sheep) (Fig. 2 A,B), summer trekking, and winter ski tourism, including artificial snow production and snow-grooming machinery that compacts the snowpack and shortens the vegetation period. Reconstruction works on the ski slopes, including the construction of water reservoirs (Suppl. material 1: fig. S1), constitute additional sources of environmental pressure.

The climate is cool-temperate, with a mean annual temperature of 5.2°C and mean annual precipitation of 1,182 mm; monthly temperatures range from −6.5°C (January) to 16.8°C (August), and precipitation is concentrated from April to August (650 mm), with May being the wettest month (166 mm) (Fick and Hijmans 2017).

The topography reflects Paleogenic folded structures. Brown forest soils dominate the montane and subalpine belts; under semi-natural pastures these degrade into meadow soils. At and above the tree line soils are diverse, including mountain meadow, mountain meadow peat, mountain meadow soddy, and mountain meadow soddy-peat (Nakhutsrishvili et al. 2006).

The vegetation zones include mixed conifer–broadleaved forests, broadly corresponding to four EUNIS (FloraVeg.EU, version 2025-10-03) habitat types: T35 – temperate continental Pinus sylvestris forest; T31 – temperate mountain Picea forest; T17 – Fagus forest on non-acid soils; and T1C – temperate and boreal mountain Betula and Populus tremula forest on mineral soils. These forests are characterised by Pinus sylvestris var. hamata Steven, Fagus orientalis Lipsky, Picea orientalis (L.) Peterm., and Betula pubescens var. litwinowii (Doluch.) Ashburner & McAll. Shrublands of Rhododendron caucasicum Pall. and Vaccinium spp. align with EUNIS S22 – alpine and subalpine ericoid heath. Subalpine and alpine pastures are analogous to EUNIS R4 ­– alpine and subalpine grasslands, with subalpine vegetation characterised by Aconitum orientale Mill., Cephalaria gigantea (Ledeb.) Bobrov, Betonica macrantha K.Koch, Lomelosia caucasica (M.Bieb.) Greuter & Burdet, and Festuca varia Haenke, and alpine vegetation typified by Festuca ovina L., Ranunculus breyninus Crantz, Alchemilla erythropoda Juz., and Sibbaldia parviflora Willd. (Nakhutsrishvili et al. 2006; Chytrý et al. 2024).

Field sampling was carried out on the Didveli ski piste, which remained accessible by cable car throughout the summer. Phytosociological plots were established within the ski run and paired with control plots in adjacent, undisturbed pasture (within 50–100 m of the piste boundary). Slope angles ranged from 10° to 45° between 1,887 and 2,664 m a.s.l. Sampling was performed during the active growing season following the methods of Braun-Blanquet (1964) and van der Maarel (2007), using the modified Braun-Blanquet scale (r, +, 1, 2m, 2a, 2b, 3, 4, 5) to estimate species cover.

The study area was divided into three altitudinal sections: (1) the highest section (Panorama piste), corresponding to the alpine zone; (2) intermediate slopes (Deka and Slalom pistes), corresponding to the subalpine zone; and (3) lower slopes (comprising Didveli 1, 2, 3, and Nodo), corresponding to upper mountainous zone (Suppl. material 1: fig. S2.). A further difference of the lower sites to the two other zones, was that all of them received artificial snow. Two parallel transects were established: one within the piste, and one immediately adjacent but outside it. Where pistes were crossed by maintenance roads, evident erosion led to the creation of a third, eroded category.

The main vegetation types within the sampled area were broadly classified according to the revised EUNIS Habitat Classification (version 2025-10-03; Chytrý et al. 2024), based on observed vegetation characteristics and floristic composition. Alpine and subalpine grasslands broadly correspond to EUNIS R4, whereas eroded and drainage-influenced areas were treated as open and disturbed vegetation types rather than assigned to a specific EUNIS habitat category.

In total, 62 plots were sampled: 26 control plots, 18 track plots (piste without additional disturbance), and 18 eroded plots (piste with additional disturbance from maintenance vehicles) (Fig. 3A-C). One control plot at 2,599 m a.s.l. was identified as a compositionally distinct outlier during preliminary exploratory analyses. This plot was strongly dominated by Rhododendron caucasicum, with only two additional species occurring at very low cover, and it was therefore excluded from subsequent analyses, leaving 61 plots in all.

Soil conditions were characterised by two variables: bare soil surface (%, assessed visually) and soil depth (cm, measured with a 40 cm auger). Plant species were identified using the multi-volume Key to Plants of Georgia and ’Flora of Georgia’ (Ketskhoveli et al. 1964, 1969; Ketskhoveli et al. 1971-2011). Specimens of uncertain identity were collected, dried, and deposited at the herbarium of the Botanical Institute of Ilia State University (TBI) for later identification. Nomenclature follows the Flora of Georgia checklist (Davlianidze et al. 2018) harmonised with international databases (GBIF 2026).

Aboveground biomass (AGB) was harvested during the summers of 2020, 2021, and 2023, at the end of July when biomass reaches its seasonal peak. Within each 25 m² phytosociological plot, a 0.1 m² frame was placed at four random positions and all above ground plant material clipped at ground level, yielding a total harvested area of 0.4m² per plot (Suppl. material 1: fig. S3). Biomass per square metre was calculated by multiplying the harvested mass by 2.5.

Harvested material was sorted into three functional groups: graminoids (Poaceae, Cyperaceae and Juncaceae), legumes (Fabaceae), and forbs (all remaining families). Material was sun-dried in sealed bags for one month at 25–38°C and weighed on a precision balance. A moisture-correction procedure was applied following Scurlock et al. (2002): a subsample from each collection was oven-dried at 70°C to constant mass, and the ratio of oven-dry to sun-dry mass was used as a correction factor to express all biomass values on an oven-dry basis (Perez-Harguindeguy et al. 2016).

Relationships among environmental variables (elevation, slope angle, bare soil surface, soil depth) and aggregate vegetation characteristics (species richness, AGB) were first explored with Principal Component Analysis (PCA) in 'between-group' mode, using Pearson's correlation coefficient and grouping plots as Control, Track, or Eroded. The resulting PCA components served as composite independent variables in a subsequent Canonical Correspondence Analysis (CCA) based on Bray–Curtis dissimilarity. PCA and CCA were performed in PAST 5 (Hammer et al. 2001). Significance of CCA axes was assessed by 999-permutation test.

Species distributions were further examined with hierarchical cluster analysis. Pairwise Jaccard distances were computed from a presence–absence matrix and plots clustered using Ward's minimum-variance linkage, implemented in PC-ORD 5.0 (McCune and Mefford 2006). The dendrogram was divided into four clusters. The relationship between cluster membership and plot categories (Control, Track, and Eroded) was tested with a Pearson chi-square test and a Monte Carlo p-value based on 2,000 replicates (χ²=21.1, p=0.0005). Clusters were interpreted based on species composition and ecological traits.

Indicator Species Analysis (ISA) was performed with the R package indicspecies (De Cáceres t al. 2010; R Core Team 2025) (multipart function, IndVal.g association measure, 999 permutations) to identify species characteristic of each plot category. Of the 211 species recorded, 156 were retained for ISA after removing singletons and invariant taxa. Species with IndVal≥0.6 and p≤0.05 were considered strong indicators; species with IndVal≥0.4 and p≤0.05 are reported as moderate indicators. These thresholds were established based on prior applications of Indicator Species Analysis (Danielsen et al. 2025; Lasram et al. 2015; Wu and Pinzon 2022) and served as pragmatic criteria rather than as a universal standard.

The impact of skiing on functional-group biomass was assessed with a further PCA (PAST 5). Year-to-year variation in functional-group biomass was tested with linear mixed-effects models (R package nlme; Pinheiro and Bates 2022), with year and land-use type as fixed effects and plot as a random effect.

Sampling yielded three plot subsets: Control (plots adjacent to but outside pistes), Track (piste plots without additional disturbance beyond skiing) and Eroded (piste plots additionally disturbed by maintenance vehicles). In total, 211 species were recorded; 103 species occurring in at least five plots were retained for multivariate analyses (Suppl. material 2: table SS1).

The between-group PCA produced two components with eigenvalues >1 (Fig. 4). PC1 (eigenvalue 4.01) explained 66.91% of variation and was strongly associated with Bare Soil Surface (negatively) and Soil Depth, Species Richness, and AGB (positively); it therefore represents a gradient of human impact intensity. PC2 (eigenvalue 1.99) explained the remaining 33.09% and was dominated by Elevation and Slope Angle, representing natural environmental variation (Table 1).

Track and Eroded plot centroids were closely associated with Bare Soil Surface, while the Control centroid was strongly associated with Soil Depth, Species Richness, and AGB. Although Track and Eroded plots overlapped considerably, both were clearly separated from Control plots, especially along PC1.

Consistent with this pattern, Control plots had substantially higher vegetation cover (mean 85%; maximum height 68 cm) and greater aboveground biomass than piste plots (Track and Eroded combined: mean cover 45%; maximum height 44 cm; minimum cover occasionally as low as 15%).

The CCA using PCA components as composite environmental variables was highly significant for both axes (p < 0.001, 999 permutations; Fig. 5; Suppl. material 1: fig. S4). CCA Axis 1 (68.62% of variation) separated Control plots from Track and Eroded plots, with the Control centroid associated with PC1 (impact-intensity proxy). CCA Axis 2 (31.38%) reflected natural environmental variation (PC2). The overall pattern of species distribution thus mirrored that of the environmental PCA: a dominant axis of human impact separating undisturbed from disturbed plots, and a secondary axis of elevational and topographic variation.

Cluster analysis identified four floristically distinct groups (Suppl. material 1: fig. S5). Chaining was low (2.02%), and cluster membership was significantly associated with plot category (Control, Track, and Eroded; χ² = 21.10, simulated p = 0.0005), indicating that the identified floristic groups were strongly related to the disturbance categories.

Cluster 1 mainly included Control plots (13 of 18 plots, 72%), representing undisturbed habitats. Cluster 2 contained predominantly Eroded (11 plots) and Track (10 plots) plots, reflecting moderate-to-high disturbance without a clear separation of erosion from track effects. Cluster 3 was compositionally heterogeneous (9 Control, 4 Eroded, 4 Track plots), corresponding to transitional communities at higher elevations where harsh abiotic conditions exert a strong compositional filter. Cluster 4 was small and predominantly Control-associated (3 Control, 1 Eroded, 1 Track). Together, the clusters confirmed that disturbance is a primary determinant of compositional differentiation, with a secondary contribution from natural environmental variation at high elevations.

ISA was performed on the three plot categories (Control, Track, Eroded), combining Track and Eroded given their strong ordination overlap. The analysis identified clear floristic differences between undisturbed and disturbed plots (Table 2).

Control plots were characterised by Gentiana septemfida (IndVal 0.683), Bromus variegatus (0.670), Veratrum lobelianum (0.636), Bistorta carnea (0.628), and Silene saxatilis (0.585). Eroded plots were associated with Draba hispida (IndVal 0.451) and Arctium lappa (0.427), whereas Track plots were characterised by Calamagrostis pseudophragmites (IndVal 0.430). Two species — Rhinanthus major var. apterus (IndVal 0.722) and Anthoxanthum odoratum (0.690) — were identified as indicators shared between Control and Track plots.

The PCA of functional-group biomasses identified two components (eigenvalues 3.7 and 1.3, explaining 74.1% and 25.9% of variation, respectively; Fig. 6; Table 3). PC1 was positively associated with Soil Depth and the biomasses of all three functional groups (especially Graminoids), indicating that deeper soils support overall higher productivity. PC2 was determined by the factor Elevation, with Legume biomass correlating negatively with Elevation, indicating that legumes grow better at lower altitudes regardless of disturbance. Graminoids and forbs were particularly sensitive to reductions in Soil Depth associated with disturbance.

Overall, graminoid biomass on Control plots was 2.5 times higher than on Track plots and 3.5 times higher than on Eroded plots. Forb biomass on Control plots was three times higher than on Track and four times higher than on Eroded plots. Mean AGB was highest on Control plots (2,526 kg ha^–1), intermediate on Track plots (1,175 kg ha^–1), and lowest on Eroded plots (833 kg ha^–1).

Forb and graminoid biomasses differed significantly among land-use types (p < 0.001) and were consistently lower on both Track and Eroded pistes than on Control plots (all pairwise p < 0.05; Fig. 7). In Control plots, forb biomass increased significantly between 2020 and 2021 (p = 0.0019), likely reflecting natural inter-annual variation. Eroded plots showed no significant temporal changes (p > 0.3), indicating persistent environmental stress and limited recovery capacity. Track plots exhibited a modest increase in graminoid biomass between 2020 and 2021 (p = 0.044) that did not persist into 2023.

Legume biomass showed a different pattern: it did not differ significantly among land-use types (p = 0.515) but increased strongly over time across all categories (p < 0.001), suggesting high resilience to disturbance and possibly benefiting from nitrogen-poor soils where competition is reduced.

Our results demonstrate that disturbance due to construction (Fig. 2. C-D) and use of ski pistes substantially alters plant community structure in the Caucasus Mountain vegetation. Control plots maintained higher species richness and aboveground biomass, representing traditional semi-natural meadow communities typical of the region's long history of pastoral land use. Eroded plots had the lowest species richness and biomass, reflecting the compounded effects of soil loss, microhabitat disruption, and reduced capacity for species establishment. Track plots showed intermediate values, with both mechanical disturbance from skiing and colonisation by opportunistic taxa.

The clear separation of Control from piste plots in ordination space, combined with substantial overlap between Track and Eroded plots, mirrors patterns described in European ski areas (Wipf et al. 2005; Klačanová et al. 2024). The overlap between Control and Track plots at high elevations is consistent with the dominance of disturbance-tolerant alpine tussock species across both disturbed and undisturbed sites at extreme elevations, as documented elsewhere in the Alps and Carpathians (Körner 2003; Klačanová et al. 2024). The overall reduction in vegetation cover from a mean of 85% on Control plots to 45% on piste plots — and the halving of maximum herbaceous layer height confirm that mechanical and chemical impacts of piste construction and maintenance exceed the natural resilience capacity of these grassland communities (Roux-Fouillet et al. 2011).

Apparently there is an interaction between elevation and impact following use as ski pistes. At higher altitudes, control and track plots exhibit greater similarity compared with lower altitudes (Fig. 4). This phenomenon may be attributed to the extensive use of artificial snowmaking at lower elevations, which consequently reduces the growing season. The outcome is a decline in species diversity and biomass.

Biomass on Control plots was 2.1 – 3.0 times higher than on disturbed plots, a ratio broadly consistent with findings from alpine ski areas in the Alps and Pyrenees (Meijer zu Schlochtern et al. 2014; Barrantes et al. 2013). The strong correlation of biomass with soil depth underlines the central role of soil integrity in sustaining mountain grassland productivity (Pintaldi et al. 2017).

The indicator species clearly distinguished undisturbed from disturbed vegetation. Control plots were dominated by characteristic alpine meadow species — Gentiana septemfida, Bromus variegatus, Bistorta carnea, Silene saxatilis, and Veratrum lobelianum — whose presence reflects intact soil conditions and absence of heavy mechanical disturbance. The abundance of Veratrum lobelianum, a species favoured by long-term overgrazing in the Caucasus (Callaway et al. 2000; Tedoradze 2022), is also consistent with the centuries-long pastoral history of the area.

Eroded plots were characterised by Draba hispida and Arctium lappa, both tolerant of, or associated with, disturbed and eroded substrates (Guna 2019). Track plots were most strongly associated with Calamagrostis pseudophragmites, which thrives in compacted soils and along drainage channels that form on flat sections of piste (Skokanová et al. 2015). The occurrence of Rhinanthus major var. apterus and Anthoxanthum odoratum as shared indicators of Control and Track plots suggests that moderate piste disturbance does not eliminate all meadow-adapted species; both are broad-habitat species of traditionally managed mountain grasslands. Rhinanthus major var. apterus is a summer annual species and is therefore less affected by mechanical stress during winter, while additional vegetation gaps created by disturbance may facilitate its seed germination (Lawniczak et al. 2011; Akhalkatsi and Tarkhnishvili 2012; Tedoradze et al. 2020).

Cluster analysis revealed a continuum of compositional change rather than a sharp boundary between disturbance states. The partial merging of Control and Track clusters at high elevations, where severe abiotic conditions override land-use effects, is characteristic of alpine ecosystems globally (Sundqvist et al. 2013). These mosaics reflect the gradual nature of ecological responses to skiing disturbance and the difficulty of establishing unambiguous reference conditions in landscapes shaped by thousands of years of pastoral management (Kikvidze 2020).

The differential responses of functional groups to disturbance provide insight into underlying mechanisms. The sharp decline of graminoid and forb biomass on disturbed plots reflects their sensitivity to soil compaction, nutrient loss, erosion, and disruption of belowground biomass by grooming machinery. Alpine forbs are characteristically slow-growing with limited regenerative capacity (Myerburg et al. 2011), and although graminoids are generally considered more tolerant of disturbance (Cole and Monz 2002), the degree of impact on ski pistes — particularly on eroded sections — clearly exceeds the resilience threshold of these communities.

The resilience of legumes is a notable finding. Legume biomass was statistically indistinguishable among land-use types and increased across all categories over the three-year monitoring period, consistent with their ability to fix atmospheric nitrogen and their competitive advantage in nutrient-poor soils (Vitousek et al. 2002; Olsen et al. 2013). Nitrogen limitation is a prominent feature of high-altitude soils (Barberis et al. 2023), and legumes are frequently among the first colonisers of degraded or disturbed alpine substrates. The elevated legume biomass on piste plots could also partly reflect deliberate sowing of legumes during piste restoration, as is practiced in some European ski areas. In contrast to our findings, Barrantes et al. (2013) reported reduced legume cover on Pyrenean ski pistes, which they attributed to very intense mechanical disturbance exceeding the tolerance even of this resilient functional group. The limiting factors for legumes may include the levels of phosphorus in the soil, as indicated by Ren et al. (2017). Furthermore, these factors — mechanical disturbance and the soil nutrient content — are interdependent.

Temporal dynamics further highlighted the limited recovery capacity of disturbed piste vegetation. The significant increase in forb biomass on Control plots between 2020 (62.4 mm) and 2021 (85.1 mm) most likely reflects natural inter-annual variability (Fig. 7). The absence of any significant change on Eroded plots over the study period indicates ongoing environmental stress. The short-lived graminoid increase on Track plots in 2021 suggests that brief episodes of partial recovery can occur but are insufficient for sustained restoration, consistent with the multi-year recovery trajectories documented in experimental ski piste restoration studies in the Alps (Wipf et al. 2005; Roux-Fouillet et al. 2011).

Grazing, which takes place on and adjacent to the ski pistes throughout summer, likely interacts synergistically with skiing impacts: piste vegetation already weakened by soil compaction and erosion is particularly vulnerable to herbivory and trampling by livestock (Mayer et al. 2009). This combination of pressures may further delay recovery and promote the dominance of disturbance-tolerant and grazing-resistant species.

This study has several limitations and methodological challenges that should be considered when interpreting the results. First, although control plots were situated 50–100 m from the piste boundary, a distance similar to or greater than that used in previous studies of ski-piste vegetation (Wipf et al. 2005; Roux-Fouillet et al. 2011), potential indirect effects extending beyond the piste boundary cannot be completely excluded. Moreover, the spatial extent of such effects may vary depending on local topography and habitat type. For example, disturbance effects may propagate differently across open grasslands than in areas where ski runs cross or border forest belts, where vegetation structure, slope configuration, and hydrological connectivity may influence the extent of impact. Future studies should therefore incorporate multiple control distances and spatial scales while accounting for topographic and habitat differences to better characterise the spatial extent of ski-piste impacts.

The temporal scope of the study represents a further limitation. Although vegetation cover and species composition were documented during all monitoring years, including 2022, biomass was not measured in that year. In addition, the relatively short monitoring period limits our ability to distinguish persistent vegetation recovery trends from interannual climatic variability and other ecological fluctuations. Consequently, the temporal patterns observed in this study should be interpreted as short-term responses rather than conclusive evidence of long-term recovery trajectories. Longer-term monitoring, ideally incorporating repeated assessments of vegetation composition, cover, and biomass, would provide a stronger basis for evaluating the persistence and stability of vegetation recovery under ski piste management.

This study provides the first assessment of ski piste impacts on vegetation in the Caucasus, filling a geographic gap in the literature on mountain tourism ecology. Our results demonstrate that a clear gradient of impact intensity exists from undisturbed adjacent pastures through standard piste tracks to heavily eroded piste sections. Ski piste construction and maintenance significantly reduced species diversity, aboveground biomass, vegetation cover, and structural complexity in montane, subalpine, and alpine vegetation, with the strongest effects on eroded piste sections. Plant functional groups responded differentially to disturbance: graminoids and forbs declined sharply, while legumes showed resilience across all conditions, likely owing to their nitrogen-fixation capacity. Grazing, occurring on and near the ski slopes throughout the summer, likely interacts synergistically with the impacts of skiing. Under the existing land-use conditions, temporal monitoring suggests relatively slow and limited vegetation recovery on disturbed piste plots during the period examined, although longer-term monitoring is needed to determine whether these patterns persist over time.

The semi-natural grasslands of the Caucasus represent a unique cultural and ecological heritage: they have been shaped by thousands of years of pastoral management (Kikvidze 2020), support exceptional plant diversity (Wiesmair et al. 2017), and form a distinctive element of the Caucasian mountain landscape (Elizbarashvili et al. 2021). The ongoing expansion of ski tourism in Georgia therefore requires careful ecological assessment and adaptive management. We recommend that future piste development minimise soil disturbance, prioritise ecological restoration of eroded sections using locally sourced seed mixtures and incorporate long-term vegetation monitoring as a standard component of ski resort management. Furthermore, it is imperative to continue with and broaden research endeavours across diverse ski slopes and resorts within Georgia to evaluate spatial and ecosystem-level variations in ecological impacts. This includes identifying vulnerable habitats and formulating management strategies tailored to local conditions. Such scholarly efforts will enhance comprehension of degradation processes and facilitate the development of effective restoration interventions. These measures are essential for reconciling the economic importance of ski tourism with the conservation of the Caucasus mountain biodiversity.