The inter-specific plant interactions along grazing and aridity stress gradients represent a major research issue in plant ecology. However, the combined effects of these two factors on plant-plant interactions have been poorly studied in the northeast of Iran. To fill this knowledge gap, 144 plots were established in 12 study sites with different grazing intensity (high vs. low) and climatic characteristics (arid vs. semi-arid) in northeastern Iran. A dominant shrub, Artemisia kopetdaghensis, was selected as the model species. Further, we studied changes in plant life strategies along the combined grazing and aridity stress gradients. In this study, we used relative interaction indices calculated for species richness, Shannon diversity and species cover to determine plant-plant interactions using linear mixed-effect models (LMM). The indicator species analysis was used to identify the indicator species for the under-canopy of shrub and for the adjacent open areas. The combined effects of grazing and aridity affected the plant-plant interactions and plant life strategies (CSR) of indicator species. Artemisia kopetdaghensis showed the highest facilitation effect under high stress conditions (high grazing, high aridity), which turned into competition under the low stress conditions (low grazing, low aridity). In the arid region, the canopy of the shrub protected ruderals, annual forbs and grasses in both high and low grazing intensity. In the semi-arid region and high grazing intensity (low aridity/high grazing), the shrubs protected mostly perennial forbs with C-strategy. Our findings highlight the importance of context-dependent shrub management to restore the vegetation damaged by the intensive grazing.

Study area

We selected two sites along the precipitation gradient in northeastern Iran in the Khorassan-Kopet Dagh floristic province of the Irano-Turanian region, located between 35°43’-36°44’N and 58°40'−60°27’ E. Based on meteorological data, Khaje has an arid climate and Baharkish has a semi-arid climate, expressed by de Martonne aridity index (see Table 1 for more details). The mean annual precipitation (20-year mean) is 255 mm in Khaje and 385 in Baharkish. A. kopetdaghensis Krasch.M.Pop. & Linecz ex Poljak was the dominant native shrub species in both sites. Artemisia species have been documented to facilitate common annual and perennial forbs (Reisner et al., 2015) by creating suitable microclimate, reducing evapotranspiration (e.g., Holthuijzen & Veblen, 2015), mediating soil temperatures (Davies et al., 2007), raising soil water content via hydraulic lift (e.g., Holthuijzen & Veblen, 2015) and accumulating soil nutrients (Cardon et al., 2013).

Sampling design

The two studied regions were 1600 ha and 1035 ha in size for the arid and semi-arid regions, respectively. The HG and LG sites were of similar size in both climatic regions. The distance between individual sampling areas within each climatic region was less than one kilometer. The HG and pairwise LG sites were relatively homogenous in terms of topography, land use, and vegetation and the only substantial difference between the paired HG and LG sites was in their grazinig intensities.The LG sites were located within fences that have prevented grazing for around 35 years, whereas HG sites were open and therefore have suffered from long-term overgrazing. Each plot was characterized by geographic coordinates and altitude. In 2017, the number of individuals and percentage cover of all vascular plant species was recorded between April and June, when the growing season peaks in this region.

The decision about the grazing status of the sites (high grazing intensity vs. occasional/low intensity grazing) was based on the median number of dung droppings (Peace 2001): 55.3 dung droppings per square meter in the HG and 6.2 in the LG sites, and also on the width of the microterrace livestock paths in a horizontal way (0.27±0.09 m for the HG site and 0.04±0.03 m for the LG site  (see more information on the grazing history in Table 1).

The sampling design was arranged in a hierarchical way: In each of the two climatic regions (arid and semi-arid), we selected six sampling areas, with a high grazed and a low grazed site in each sampling area, arranged in a pairwise way (hereafter referred to as HG and LG sites). Then, we sampled three plots under the A. kopetdaghensis shrubs and three adjacent plots outside the canopy of A. kopetdaghensis (hereafter referred to as under-canopy and open plots) in each HG as well as LG site (Soliveres et al. 2014b). Altogether, 144 plots were sampled: 2 climatic regions, 6 sampling areas in each climatic region, one pairwise HG and one LG site in each sampling area and 6 plots (3 under-canopy and 3 in the open) in each HG or LG sites (see Figure 2). We recorded the numbers of individuals of all vascular plant species and their percentage covers and then calculated the Shannon index of species diversity (H = -∑ pi ln pi) for each plot (Shannon, 1948); pi is the proportion (n/N) of individuals of one particular species (n) divided by the total number of individuals (N).

To obtain comparable samples in the surrounding ‘open’ plots (outside the canopy of A. kopetdaghensis), matching the size of each sampled A. kopetdaghensis canopy, we sampled at randomly selected paired points, located  ~1 m away from the canopy edge of each sampled A. kopetdaghensis shrub. When the size of A. kopetdaghensis was not measured, a wire loop was shaped to match the size of the sampled A. kopetdaghensis canopy plot and then used to define the size of the patch sampled in the ‘open’ plot (Farzam & Ejtehadi, 2017). In addition, percentage covers of all vascular plants in plots in these open areas were recorded and identified to the species-level.

Statistical analyses

Relative interaction intensity (RII) was used to assess the effect of shrubs on under-canopy vegetation (Armas et al., 2004) and was calculated based on the cover, richness, and diversity (expressed as Shannon index) of under-canopy vegetation: RII = (value under shrub – value in the open)/(value under shrub + value open). Samples were paired between each A. kopetdaghensis shrub and its neighbouring open plot. RII was used as an indicator of the facilitation by the target shrub, based on the performance of under-canopy plants. The interaction index has defined limits [-1,+1], with positive values indicating facilitation and negative values indicating competition.

The differences in RII indices for species richness, cover, and diversity between the HG and LG sites and between the arid and semi-arid regions were tested using linear mixed-effect models, with “sampling areas” as a random effect, “climatic region” and “grazing” as fixed effects and RII based richness (RII-Richness), cover (RII-Cover), and Shannon H (RII-Shannon diversity) as response variables. All univariate analyses were performed in the R software, using the NLME package. The script for the model testing the interaction between “climate” and “grazing” were “lme(Relative interaction intensity~climatic region*grazing, random=~1|sampling area)”. The normality of the input data was assessed based on Shapiro-Wilk tests, and the normality of residuals was checked visually, by plotting the observed values against the fitted values.

Further, we used the method of indicator species analysis to reveal the preference of individual species for the HG versus LG sites in both the arid and semi-arid climatic regions. With this approach, we could determine the indicator species sensitive or resistant to high grazing intensity in two different climatic regions. Indicator species analysis has two main components: (i) recorded on either HG or LG sites only (exclusivity); (ii) recorded on all samples of either the HG or LG group (fidelity). The indicator value index was assigned to all species, identifying species with the highest association values. The permutation tests (999 permutations) were used to estimate the statistical significance of individual species’ indicator values (Dufrêne & Legendre, 1997). The indicator species analyses were performed using the “indicspecies” package of the R software (R Development Core Team, 2013).

We also calculated the values for CSR plant strategies for all indicator species as well as for A. kopetdaghensis, following Pierce et al. (2017), based on the following traits: specific leaf area (SLA), leaf dry matter content (LDMC) and leaf area (LA). We collected the leaves from robust and well-grown plants. Leaf material was collected from 10 individuals of each species (Behroozian et al., 2020), packed in moist paper bags, sealed in plastic bags, and stored in a thermal box until storage at 4 °C for 12–24 h. Depending on the size of leaves, 2–10 undamaged, fully expanded young leaves (including the petiole) were measured per individual. We determined the leaf area using a digital scanner and Leaf Area Measurement v1.3 software (Andrew Askew, University of Sheffield, UK). Turgid leaf fresh weight (LFW) was obtained from saturated leaves, and leaf dry weight was determined after drying for 72 h in an oven at 70 °C. For CSR strategy analysis, values of LA, SLA, and LDMC were inserted into the ‘StrateFy’ spreadsheet 3 to calculate C, S, and R percentages for each species (Pierce et al., 2017).

The datasheet including the climate (Arid/semi-arid), sampling area, grazing (grazed/ungrazed), and plant-plant interaction (the plot under shrub canopy/open area) treatments. Also, the plots contain the present % cover of each species.