Changes in Semi-Arid Plant Species Associations along a Livestock Grazing Gradient

In semi-arid ecosystems, vegetation is heterogeneously distributed, with plant species often associating in patches. These associations between species are not constant, but depend on the particular response of each species to environmental factors. Here, we investigated how plant species associations change in response to livestock grazing in a semi-arid ecosystem, Cabo de Gata-Níjar Natural Park in South East Spain. We established linear point-intercept transects at four sites with different grazing intensity, and recorded all species at each point. We investigated plant associations by comparing the number of times that each pair of species occurred at the same spatial point (co-occurrences), with the expected number of times based on species abundances. We also assessed associations for each shrub and grass species by considering all their pairs of associations and for the whole plant community by considering all pairs of associations on each site. At all sites, the plant community had a negative pattern of association, with fewer co-occurrences than expected. Negative association in the plant community increased at maximum grazing intensity. Most species associated as expected, particularly grass species, and positive associations were most important at intermediate grazing intensities. No species changed its type of association along the grazing gradient. We conclude that in the present plant community, grazing-resistant species compete among themselves and segregate in space. Some shrub species act as refuges for grazing-sensitive species that benefit from being spatially associated with shrub species, particularly at intermediate grazing intensities where positive associations were highest. At high grazing intensity, these shrubs can no longer persist and positive associations decrease due to the disappearance of refuges. Spatial associations between plant species and their response to grazing help identify the factors that organize plant communities, and may contribute to improving management of semi-arid ecosystems.


Introduction
Plant species associations are a fundamental aspect of plant community ecology [1][2][3]. Analyses of plant species associations provide information about environmental heterogeneity, biotic interactions and patterns of seed dispersal [4][5][6][7]. This information is of particular interest in semi-arid plant communities where vegetation often occurs in patches. Usually, vegetation patches are composed of shrubs that can act as shelter against harsh environmental conditions. These shrubs are called 'nurse plants', for they appear to provide microhabitats that enhance survival for other plant species [2,8]. Theoretical models based on empirical studies suggest that these positive interactions between plant species are one of the main drivers for the formation of these patches [8][9][10].
Since the article by Callaway and colleagues [11], additional studies have proliferated showing the presence and importance of positive interactions in plant communities [6,7,10]. It has been suggested that positive interactions should be particularly common in communities developing under high stress conditions and in those exposed to high consumer pressure [2,11]. This theory has been referred to as the 'Stress Gradient Hypothesis' (SGH) and can explain some of the patterns in plant species interactions occurring in stressed ecosystems (but see [12,13]). However, most studies have focused on a limited number of species within a community [14], and analyses at the community level have not provided unequivocal support for SGH [15][16][17]. Experiments at this scale are complicated because each species responds in a particular way to each stress, and, typically, those responses can change depending on the species' ontogeny, habitat, and type of stress considered [17][18][19]. Usually, the response of a community is viewed as the net importance of positive and negative interactions in the community (sensu, the proportion of the total interactions in the community that are positive or negative [20]).
One alternative to studying interactions at the community level has been using spatial association between plant species as a surrogate for quantifying interactions [7,15,16]. This correlative approach assumes spatially associated species result from positive interactions, and species that are negatively associated are segregated by competitive interactions [7]. When interactions between species are weak, plant species will associate at random. Spatial association has been employed for studying plant species interactions in arid communities. For example, Verdú and Valiente-Banuet studied positive spatial associations between shrubs and seedlings in the Sonoran Desert, and found that there was a relationship between positive interactions and co-evolutionary processes in that plant community [21].
Grazing is one of the most important biotic factors shaping plant communities. Biomass consumption by herbivores affects both plant species composition and community spatial structure [22][23][24][25]. In arid and semi-arid ecosystems, grazing reduces total plant cover, increases abundance of certain life forms such as annual plants, and changes the identity of dominant species [24,26]. Moreover, grazing may lead to increased positive interactions between plants as a result of associational defense; for example, some plant species protect themselves from herbivores by spatially associating with unpalatable plants [27][28]. Previous works testing SGH in ecosystems under grazing stress have found an increase in the importance of positive interactions at low grazing levels, but not at high grazing levels where negative interactions predominate [13,29]. Assessing the effects of grazing on plant interactions provides valuable information for ecosystem management (e.g. which species act as refuge for grazing-sensitive species; which species need a refuge to survive). Changes in community structure are central to detecting when an ecosystem is overgrazed [24].
In the present work, we analyze the spatial associations among all plant species in a semi-arid community occurring along a gradient of livestock grazing. We evaluate associations for each individual species, between all pairs of species, and in the whole plant community, and how these associations change due to grazing. We estimate all associations by comparing real spatial cooccurrences between plant species with expected co-occurrences due to species abundances. Specifically, we hypothesize that the whole plant community will become more positively associated at intermediate grazing intensities, and that associations of each plant species and between pairs of species will depend on species life forms. We distinguish among associations between shrubs (those species responsible for patch formation), between grasses, and between shrubs and grasses.

Study Area and Data Collection
The study was conducted in Cabo de Gata-Níjar Natural Park, which lies along the Mediterranean coast in Southeastern Spain (36u 469 N, 2u 099 W). The park occupies 37,570-ha park, with a maximum elevation of 493 m (El Fraile Peak). The climate is semiarid Mediterranean (marked seasonality, drought in summer and most rainfall in spring and autumn. Average annual rainfall = 193.9 mm, Average annual temperature = 19.4 C [30]). Historically, the area has been used as an agro-pastoral system, with cereal cropping on floodplains and livestock (sheep and goat) grazing on slopes all the year. The plant community is characterized by Chamaerops humilis L., and other common species include Rhamnus lycioides L., Pistacia lentiscus L., Periploca laevigata Aiton, and Stipa tenacissima L. [31].
Vegetation data were collected from the Southern section of the park, a volcanic area where highly stony soils predominate [32]. All permits required to carry out the field studies were obtained from the Natural Park authorities. In that region, vegetation is characterized by open shrubland with shrubs organized in patches, which are embedded in a matrix of a large tussock grass, S. tenacissima L. S. tenacissima is a very strong competitor that colonizes the gaps within patches caused by livestock and aridity, and can exclude other species once it becomes the dominant species [10,33].
Within the area, four sites at different distances from El Romeral farm were selected. Movements of animals (sheep and goats) were monitored for one week per season with a GPS device. Effective stocking rate was calculated from the average stocking rate of the farm (ind?ha 21 ) multiplied by a correction factor based on the percentage of time each grazing site was used. Sites were ranked based on the amount of grazing pressure to which they had been exposed (G1 = 0 ind?ha 21 ; G2 = 0.27 ind?ha 21 ; G3 = 0.46 ind?ha 21 ; and G4 = 0.65 ind?ha 21 ). Grazing carrying capacity for this plant community is 0.39-0.57 ind?ha 21 , so we considered G2 as low grazing intensity, G3 as intermediate grazing intensity and G4 as high grazing intensity [34]. In April, 2001, three 500-m-long linear transects were established at each site, and the Point-Intercept Method [35] was used to quantify vegetation. On each transect, the species occurring at each point were recorded every 20 cm. Presence and life form (shrub or grass) of all species were recorded and no distinction was made between the ontogenetic stages of individuals. All transects were run parallel to the slope, separated by 50 m and established at the same altitude, orientation and soil parent material.

Association Measurements
A plant-plant association matrix A SxS was built for each site, based on the data pooled from the three transects. S is the number of species present at the site. In the matrix, a ij is the number of times that species i and j co-occurred at the same sampling point (with a ij = a ji ). The matrix was used to calculate the total number of co-occurrences for a given species i (a i~P S j~1 a ij ) and the total number of co-occurrences at a site (A~P S i~1 a i =2). The diagonal terms of the matrix were set to 0 because it was not possible to estimate the co-occurrence of a species with itself from the presence data.
To test the deviation from the expected patterns of plant species associations, an E SxS matrix was calculated for each site. This matrix includes the expected number of co-occurrences between species based on their abundances. For each species i, its relative cover was calculated as p i = n i *T, where n i is the number of points where species i occurred at each site, and T is the total number of points sampled (2501 sampling points * 3 transects = 7503). In that context, p i is the probability of finding the species i at one randomly selected point at a site. Using those data for all of the plant species, a P SxS matrix was computed, where p ij = pi*pj. Therefore, p ij is the probability of finding the species i and j at the same sampling point on the site. The expected co-occurrences matrix E SxS was computed as E SxS = P SxS *T (with e ij = p ij *T), and total expected co-occurrences for each species i (e i~P S j~1 e ij ) and each site (E~P S j~1 e ij =2) were calculated similarly as for a i and A.
A Poisson distribution was employed to compare A SxS and E SxS . The Poisson distribution is a statistical distribution widely used for analyzing count data. Furthermore, it has long been used in ecological analyses [36], particularly with vegetation data collected using the Point-Intercept Method (i.e. the number of contacts of a given species fits a Poisson distribution if individuals are distributed randomly and the probability of more than one contact for the same individual is negligible [35]). The Poisson distribution is characterized by one parameter, l, which determines the mean and variance of the distribution. Thus, each value of A, a i and a ij was compared to a Poisson distribution fitted with its corresponding value, E, e i and e ij , as the l parameter. To determine whether observed co-occurrences (A, a i and a ij ) differed significantly from the co-occurrences expected based on species abundances (E, e i and e ij ), 95% confidence intervals were calculated for each distribution. The comparison between A and E is the general association pattern present in the plant community (plant community presents more, less or the expected total number of cooccurrences). The comparison between a i and e i is the general type of association of species i (species i presents more, less or the expected total co-occurrences with all other species), while the comparison between a ij and e ij is the particular association between species i and j (species i co-occurs more, less or as expected with species j). Plant species were believed to be positively associated when they co-occurred at a level greater than that expected by chance. Negative associations between plants were inferred when co-occurrences were less frequent than expected by chance'. If actual co-occurrence values fell within the confidence interval, these values did not differ significantly from those expected because of species abundance and were considered a neutral association.
To assess possible changes in community associations along the grazing gradient, the importance of positive/negative associations at each site was measured as the relative increase/decrease in total co-occurrences with respect to expectations (R = (A -E)/E). The positive association dominates the plant community if R .0, whereas if R ,0, the negative association is more dominant. If R , 0, neither association dominates. Total co-occurrences (A and E) include all species in the plant community, but do not provide information about the number of species or pairs of species exhibiting each particular type of association. Therefore, we calculated the proportion of species and pairs of species that presented positive, negative and neutral associations. However, when species are uncommon it is not possible to detect a negative association because 0, the minimum observable value for a i and a ij , falls within the 95% confidence interval of the expected distribution. Therefore, we only considered species and pairs of species that were sufficiently abundant so that we could distinguish between positive, negative and neutral associations. For a i , the proportions of species considered at each site were G1 = 70%, G2 = 69%, G3 = 57%, and G4 = 64%; and for a ij the proportions of pairs of species considered were G1 = 3.3%, G2 = 4.1%, G3 = 3.3%, and G4 = 3%. Because it was possible to distinguish between neutral and positive associations among all species and pairs of species, we calculated the importance of positively associated species and pairs of species as the proportion of positive associations in relation to all possible associations (species, Rs = s+/S; and pairs of species, Rss = ss+/S(S -1); where s+ is the number of species showing a positive association, and ss+ is the number of pairs of species that are positively associated with each other). As Rs and Rss increase, more plant species represent a positive association, and more pairs of species are positively associated. Rs was calculated for both life forms and for all species, and Rss was calculated for associations between species with the same life form, between different life forms and between all species. All analyses were performed using R (http://www.Rproject.org). The variables and parameters used in the analyses are presented in table 1.

Results
Grazing modified the plant community in Cabo de Gata-Níjar Natural Park. Species richness was 70% greater at the ungrazed than grazed sites, particularly due to the large number of grass species, and biodiversity decreased as grazing became more intense ( Table 2). The mean number of species recorded per point was largest at G1 and the number of points with no records (Bare soil) was largest at low and intermediate grazing values (G2 and G3). Mean number of co-occurrences per point decreased with grazing. Abundance (n i ) and co-occurrences (a i ) of each species for each grazing intensity are included as support information (Table  S1).
At all four sampling sites, plant communities exhibited fewer total co-occurrences (A) than were expected by chance (E), which indicated that plants were more likely to be alone in this community, rather than in association (Fig. 1). Negative association was most important at the highest grazing intensity (R G1 = 20.233, R G2 = 20.262, R G3 = 20.269, R G4 = 20.476).
Regarding the general type of association of species (a i ), at all sampling sites there were species which associated positively,   neutrally or negatively, but species that associated neutrally were the most common. Besides neutrally associated species, those negatively associated were more frequent than positively associated ones (Fig. 2a). No species presented a shift in its type of association along the grazing gradient (i.e. no species exhibited both a positive and a negative type of association along the gradient, Table 3). Negatively associated shrubs were common, particularly at the non-grazed site (Fig. 2b), while negatively associated grasses increased at sites with highest grazing (Fig. 2c).
For associations between pairs of species (a ij ), all three types of association were also detected. Neutral associations were the most common, and negative associations were more common than positive ones (Fig. 3a). Negative associations between shrubs were highest in areas with low and intermediate grazing intensities (Fig. 3b), whereas negative associations between grasses and between shrubs and grasses were highest at high grazing intensity ( Fig. 3c-d). Individually, each species could associate positively with some species, and neutrally or negatively with the rest. The importance of species with positive association (Rs) was highest at low grazing levels for both life forms and for all species (Fig. 4a). The importance of positive associations between all pairs of species (Rss) remained nearly constant, but they decreased at the site with the highest level of grazing (Fig. 4b). The importance of positive associations between shrubs was highest at low grazing intensities, whereas between grasses and between shrubs and grasses it was highest at the non-grazed site.

Discussion
In Cabo de Gata-Níjar Natural Park, plant community responded to grazing intensity. Grazing reduced community biodiversity, and increased bare soil except at the site with highest grazing intensity ( Table 2). Grazers preferentially feed on palatable species, favoring the persistence of non-palatable species [24]. This selective removal modifies the abundance of plant species, with some species becoming dominant while others disappear from the community [26]. The effect of grazing was particularly remarkable for grasses, which lost half their species. When grazing became very intense, the plant community was dominated by one very abundant species (S. tenacissima). These changes in species composition and abun-dance also modified the associations between species in the plant community.
Regarding the general association pattern in the community, there were fewer co-occurrences than expected by chance. At all sites, irrespective of grazing intensity, plants showed a tendency to segregate rather than associate. Thus, it seemed that negative associations dominated the community, particularly at maximum grazing intensity. In recent works in the Spanish Mediterranean region, positive interactions between plant species have been presented as a main determinant of the plant community [37]. Here the predominant interaction across all grazing intensities was negative. Furthermore, associations among abundant species were predominantly negative. Grasses such as S. tenacissima, and shrubs like Launaea lanifera and Thymus hyemalis generally associated negatively with other species independently of grazing level. These species are common in degraded habitats and, through competitive exclusion (S. tenacissima [33]) or allelopathy (T. hyemalis [38]), they usually occur alone in the area rather than in association. In this area, abundant plant species are adapted to the harsh semi-arid environment and grazing, and they compete with other well-adapted species for space and other resources [39]. On the other hand, some abundant species benefit from the association with other plant species. For example, the grass Brachypodium retussum preferentially develops under the canopy of other species [33] and the shrub Phlomis purpurea presents lowdensity branching, which allows it to enter vegetation patches [40]. Although some plant taxa exhibited positive associational patterns, neutral and negative associational patterns were most common across all species in the communities.
Large shrubs such as Chamaerops humilis, Genista ramosissima and Periploca laevigata are responsible for the development of vegetation patches in the area [31]. Often, in semi-arid environments, these shrub species act as 'nurse plants' because they facilitate establishment and development of species under their canopy [21,41,42]. In our study, several positive associations between shrubs or between shrubs and grasses reflect this nursery effect. However, negative associations between these shrubs and the competitive abundant species overcame the positive association that shrubs established with other species. There are some examples of the facilitative effect of grasses on the establishment of other species in semi-arid ecosystems [43], but in our case most of the positive associations included at least one shrub, suggesting the role shrubs have as 'nurse plants' in the community.
In order to analyze associations of species and pairs of species, we considered only those species abundant enough to allow distinguishing negative from neutral associations. Typically, in plant community studies uncommon species are excluded from the analyses because they do not provide robust results [7,44]. In our case, despite the low likelihood of detecting negative associations in uncommon species, we found many positive associations between these species. Others have suggested that rare species are more likely to be facilitated than abundant ones [44]. Our results also suggest that rare species are likely to be associated with other species in semi-arid plant communities.
In our study, the proportion of positive species associations increased at low and intermediate grazing intensities (G2-3) and was lowest at the highest intensity (G4). This result has been reported in other studies dealing with changes in the interactions between plant species along grazing gradients, but to our knowledge this is the first time that this result is evaluated at community level [13,26,29]. As grazing increases, associations between plant species become more frequent, particularly those associations with shrub species that act as refuges against grazers. However, once a particular threshold is reached, grazing intensity Chamaerops humilis* L.  is such that facilitative species cannot provide protection and positive associations decrease [29]. Interestingly, the particular type of association of each species remained consistent across the grazing gradient. One possible explanation is that each species exhibits a predominant associative type in the community, regardless of grazing intensity (i.e., competitive species never become associative [25]). Others have documented associational changes across different stages of plant ontogeny [19], but here the observed associations reflect the general type of association of the species in the area. Spatial associations between species have been presented as an indirect measure of species interactions, but this approach has limitations. The spatial association of species is the net result of biotic interactions, seed dispersal and environmental heterogeneity [45]. The present study has additional limitations. For example, inter-specific associations are not measured, all individuals of each species are considered ecologically identical (e.g. different life stages may interact differently [19,46]) and other effects of interactions are ignored (e.g changes in species biomass [14]). Nevertheless, our results are consistent with those reported by works about SGH including a limited number of species [13,27]. Identifying the processes that drive the organization of natural communities under grazing and the role that each species plays in said organization provides valuable information about which species maintain the structure of grazed ecosystems. This information is central to detect overgrazing events in grazed ecosystems [24].
In conclusion, most species were not associated with other species and the most common association among plants in this semi-arid plant community was negative, especially associations with dominant species. This suggests that either neutral processes and/or competitive interactions are structuring these plant communities. The associational patterns of most species did not vary with grazing intensity; however, there was a tendency for positive associations among species to become less frequent at high levels of grazing. Positive associations among plants appeared to be most important at low and intermediate grazing intensities. Identifying non-neutral species associations provide information about the processes and species driving the organization of natural communities and helps further the development of conservation and restoration plans.

Supporting Information
Table S1 Life form, abundance (n i ), co-occurrences (a i ) and association type (ass) of plant species in Cabo de Gata-Níjar Natural Park along grazing gradient. Association values are presented for species that could distinguish between neutral and negative associations. (XLSX)