Study site

El Potrero Chico (Nuevo León, Mexico) is located on the northern edge of the ‘Sierra el Fraile y San Miguel’ Natural Protected Area, which has an area of 23,506 ha between 800 and 2,360 m a.s.l. This area is part of the Sierra Madre Oriental mountain range. Rock types include sedimentary rocks of marine origin dating back Mesozoic era, shale, and limestone cliffs (INECC 2017). Limestone supports abundant and diverse cliff vegetation, and favoured the evolution of plant specialization (Farris 1995; Larson et al. 2000; Lichter-Marck and Baldwin 2023). El Potrero Chico has a semi-arid climate with hot summers (average monthly maximum temperatures over 35 °C between June and August) and moderate cold temperatures during winter (average monthly minimum temperatures between 7 and 16 °C). The highest precipitation is reached in September and October with averages ranging from 70 to 130 mm while the rest of the year monthly precipitation is below 50 mm.

El Potrero Chico is one of the world’s prime climbing destinations with over 600 climbing routes grouped into 24 climbing sectors. The first recorded climbing in El Potrero Chico was in 1960 but climbing experienced rapid development in this area from the late 1980s. Winter and early spring (between November and May) are the seasons when most climbers visit the area due to the moderate temperatures and lower precipitation. The number of climbers is lower from June to October due to the hot or rainy weather.

Sampling design

To examine differences between unclimbed and climbed routes, we designed a closely-adjacent case-control sampling design with a 3 m wide × 3 m high quadrat placed along the climbing route (Fig. S1). The quadrat was composed by a central Climbed (C) plot of 1 m wide and 3 m high, two immediately adjacent plots of 0.5 m wide and 3 m high, which were not surveyed (i.e. no data were obtained from these plots), and two Unclimbed (U) plots of 0.5 m wide and 3 m high on the left and right side of the 3 m × 3 m quadrat that were used as controls, since they represent areas not reached by climbers (Fig. S1). The use of a closely adjacent paired design is essential to adequately test the impact of rock climbing on cliff vegetation (Boggess et al. 2021), since this precludes the possibility that variations in biotic or abiotic factors such as aspect, micro-topography and insolation that could act as drivers of differences between climbed and unclimbed plots (Holzschuh 2016; Boggess et al. 2021). Closely-paired transects have the added benefit of avoiding an observer’s interference in the undisturbed areas, since unclimbed transects can be surveyed from the same anchor with the help of directional gear placements (Boggess et al. 2021).

To define the position of the climbed plots (and thus, of the sampling quadrat), the bolts installed in the cliff-face were considered as the central point (i.e. 0.5 m to the right and 0.5 m to the left of the bolt), since the bolt represents with high precision the typical middle point that climbers use when ascending. However, to avoid interference with adjacent climbing routes, the selected routes for sampling were at least 5 m distant from the next climbing route. The unsurveyed plots are adapted from March-Salas et al. (2018) and guarantee separation between the unclimbed and climbed plots, since not all climbers follow exactly the same path along a climbing route. This prevents biased data acquisition from casual climber’s ascent deviations, as unsurveyed areas cannot be considered completely undisturbed (Boggess et al. 2021). Moreover, in order to characterize the spatial distribution of plants within each plot, both climbed and unclimbed plots were divided in 0.5 m × 0.5 m subplots (i.e. 12 subplots in each climbed plot and 12 subplots in each unclimbed plot; see Fig. S2). Pictures were taken from each subplot (see below in ‘Data collection’).

To examine the maximum spatial distribution of cliff-face plants, we established the sampling quadrats at three heights along the climbing route, positioned at the Top, Middle and Bottom (Fig. S1). In order to fit the three quadrats without any overlap, we selected climbing routes of between 15 and 35 m height. The distance from the Middle to the Top and Bottom plot was roughly equidistant.

Data collection

Field surveys were conducted from November 2019 to December 2020. We sampled 12 climbing routes of El Potrero Chico (Table S1), adding up 36 climbed plots with a sampled area of 108 m2 and 36 unclimbed plots with a sampled area of 108 m2. The sampled routes were placed in contrasted aspects: North (n = 5), South (3), East (1) and West (3). We noted the height of each climbing route as well as the climbing difficulty using the Yosemite Decimal System (YDS), grouped in three classes in our sampling sites: beginner (5.6–5.9), intermediate (510a–5.11d), advanced (5.12a–5.13d). To account for physical characteristics of the rock (i.e. micro-topography), we measured the slope of the center of each quadrat in the field, and the proportion of cracks (i.e. crevices) in each 0.5 m × 0.5 m subplot using ImageJ, and the estimated both measures at the plot level. These measurements are crucial to eliminate potential bias when testing the climbing effect, since the establishment and survival of plants is more restricted under steeper and negative slopes, and with a lower abundance of cracks (Larson et al. 2000; Holzschuh 2016).

To determine the climbing intensity of each route, we used the Climbing-Use Intensity (CUI) index developed by Clark and Hessl (2015), as a function of the walking time required to reach the cliff base and the popularity of the climbing route inferred by the number of stars (0–4) assigned in a reference and updated climbing guidebook of the area (Madden 2019). In order to use a standardized and categorized measure, we grouped the CUI values by quartiles (Clark and Hessl 2015), resulting in low (Q1), moderate (Q2), high (Q3) or very high (Q4) climbing intensity.

We noted all the species of vacular plants to calculate the species richness in the climbed and unclimbed plots of each route and quadrat, as well as the number of individuals per species (i.e. abundance). Unidentified species in the field were later identified through image determination by local botanical experts but 18 of the 63 species could only be determined at the genus level. Species were further classified as endemic (species restricted to the Sierra Madre Oriental, Mexico), native (non-endemic but present in Mexico), and alien species (Velazco et al. 2011; Salinas-Rodriguez et al. 2017), and according to their rock association as rock-specialists (i.e. restricted to rocky habitats), species with non-strict but close association to rocky habitats (i.e. frequently inhabiting rocky environments but also found in other ecosystems) and generalist species (see Table S2). Shannon-Wiener (H’) and Simpson (D) diversity indices were calculated per cliff as well as for climbed and unclimbed plots within each route using the diversity function from the vegan R package (Oksanen et al. 2020). Based on photos, plant cover was determined by the area (i.e. plant orthogonal projection) using ImageJ (in cm2). We then calculated the percentage of a plant’s cover relative to the size of the climbed or unclimbed plot. Additionally, the relative cover (CR i ), the relative abundance (AR i ) and the relative frequency (FR i ) of each species in the sampled plots were calculated (Alanís et al. 2020). We also classified species into dominant (DO), common (CO) and locally rare (RA) by using the Importance Value Index (IVI) of species together with the species distribution range and local presence (Curtis and McIntosh 1951; Velazco et al. 2011). IVI was calculated using the importance value function from the Biodiversity R package (Kindt and Coe 2005). IVI considers the sum of the relative frequency (number of plots where a species is observed divided by the total number of surveyed plots), the relative abundance (in terms of number of individuals of a species, also referred to as relative density) and the relative spatial dominance (in terms of percentage of rock area cover by a given species) of species. These calculations determine the ecological value in terms of abundance and biomass and thus the dominance of the species in the plant’s community (Curtis and McIntosh 1951). Species with the 15% highest IVI were considered dominant species (DO), species with the 15% lowest IVI were considered locally rare species (RA), and species with in-between IVI values were considered common species (CO, Table S2).

Data analysis

We conducted all statistical analysis with R version 4.0.3 (R Development Core Team 2020). We used Linear Mixed-effects Models (LMMs) implemented in the lme4 package and the lmer function (Bates et al. 2015) to test the effect of rock climbing (referred to as climbing) on plant abundance, cover, and species richness, and whether this effect differed among different climbing intensity levels. Plant abundance, cover, and species richness were included as response variables in three separate models. Cliff section (three levels: Bottom, Middle, Top), climbing difficulty (three levels: beginner, intermediate, advanced), climbing effect (two levels: climbed vs. unclimbed), climbing intensity (four levels: low, moderate, high, very high) and the two-way interaction between climbing effect and climbing intensity were modelled as fixed factors. Climbing route nested in climbing sector (i.e. a climbing area with multiple routes) was included as random factor, and the slope and the percentage of cracks as covariates. Additionally, in two separate models, we used LMMs including climbing effect as fixed factor and route nested in sector as random factor to test whether Shannon-Wiener and Simpson diversity indices calculated per study site (i.e. route) differed between climbed and unclimbed plots.

To detect patterns of co-occurrence among cliff species and whether this co-occurrence varies among the species-dominance levels (i.e. rare, common, dominant species), we used the cooccur function from the cooccur R package (Griffith et al. 2016). This species co-occurrence analysis was conducted for species occuring in the same route, same cliff section and same climbing effect, as interaction would occur at this spatial level. In this way, the presence/absence co-occurrence matrix of all species (Fig. S3) and the co-occurrence within climbed and unclimbed plots were analyzed and mapped in different figures, highlighting the positive, negative (both considered as non-random associations) or random associations. If the presence of one species favours the presence of another in a non-random way, the association is considered positive; the association would be negative if the presence of one species systematically hinders the presence of another. Random associations are those that do not deviate from their expected co-occurrences by more than 0.1 considering the number of plots generated (Griffith et al. 2016). Only co-occurring species are shown in the matrix, so the analysis represents an approach of the number of species co-occurring, and thus coexisting and potentially interacting in each condition (i.e. by climbing effect and route section). Subsequently, we calculated the number of co-occurrences between pairs of groups of species dominance level (as explained above, classified according to the Importance Value Index – IVI – of the sampled species; Curtis and McIntosh 1951; Velazco et al. 2011). We also used LMMs to test how climbing affected the abundance and cover of each of the three dominance groups. These models included the group of species dominance level (three levels: rare, common, dominant), climbing effect and their two-way interaction as fixed factors, and route nested in sector as a random factor.

Finally, we tested for changes in community composition between sites that can be attributed to climbing. To this aim, we first used permutational multivariate analysis of variance using distance matrices with the adonis function from the vegan R package (Oksanen et al. 2020) in order to assess the extent that factors influence the species composition while controlling permutations by routes (i.e. sites). Second, we implemented non-metric multidimensional scaling (NMDS) analysis to visualize and interpret the species configuration according to climbing. We also conducted NMDS analysis for testing for variation in species composition among and within communities (i.e. among routes, and among the bottom-middle-top sections within each route). Here, we used the MetaMDS function of the vegan R package (Oksanen et al. 2020) that calculates Bray-Curtis distances for the community-by-site matrix.

In all LMMs, we tested the assumptions of normality and homogeneity of variance of the residuals using the Shapiro-Wilk test and the Bartlett test, respectively, and also checking visually. If the residuals were not normally distributed, we transformed the response variable. In the case of heteroscedasticity, we applied a weighted least square regression (Strutz 2016) by including weights (1/variance) into the model, using the extract model weights command. Whenever there were significant main effects containing more than two levels or significant interactions, we applied post-hoc contrasts using the lsmeans package (Lenth 2016) with the Tukey’s test.