Increased grazing drives homogenisation but reduced grazing increases turnover in upland habitat mosaics

Many ecosystems are grazed by livestock or large, wild herbivores and exist as mosaics of different vegetation communities. Changing grazing could have an impact on heterogeneity as well as on composition. A long-term, large-scale grazing experiment that maintained existing low-intensity sheep grazing, tripled it, removed it and partially substituted sheep grazing by cattle grazing was set up on a mosaic of upland vegetation types. The impact of changing grazing regimes was assessed in terms of changes in temporal and spatial species and functional beta diversity. Removal of grazing had the highest impact on species replacement, whilst increased grazing was closest to maintaining the original species complement. Wet heath and Molina mire had the lowest turnover, but wet heath showed the highest changes in unidirectional abundance as it contained species capable of increasing in abundance in response to changing grazing intensity. Agrostis-Festuca and Nardus grasslands displayed the highest level of balanced species replacement reflecting their more dynamic vegetation. In functional terms, there was no clear separation of communities based on their grazing preference, all were relatively resistant to change but Nardus grassland was the most resistant to the removal of grazing. The increased offtake associated with increased grazing led to a degree of homogenisation as grazing tolerant species associated with preferred communities increased in the unpreferred ones. Decisions about grazing management of the uplands involve many trade-offs, and this study identified potential trade-offs between stability and homogenisation to add to existing ones on the biodiversity of different groups of species and on ecosystem services.


Introduction
Many ecosystems are grazed by domestic livestock or large, wild herbivores. Changing grazing management or herbivore population sizes can drive changes within plant communities leading to shifts in plant species abundances and diversity (Hester et al. 2006) and often to shifts in ecosystem function driven through changes in the functional traits of the species present (Milchunas et al. 1998). Many communities exist as mosaics driven by interactions between grazing and underlying differences in edaphic or topographic conditions and dominant species in one community may be subordinate species in other communities and vice versa (Rodwell 1991(Rodwell , 1992. Depending on how these dominant and subordinate species react to changes in grazing, a change can lead to increased heterogeneity or homogeneity in species composition and functional traits as abundances and distributions change in response. Changes in grazing could have an impact on heterogeneity (beta-diversity, i.e., the differentiation in species between communities, Whittaker 1960) within ecosystems. Increased grazing should favour tolerant or resistant species whilst a reduction in grazing could promote taller, grazing intolerant species and lead to the exclusion of low-growing species (Augustine and McNaughton 1998). Within a mosaic of communities that share species, increased grazing could reduce heterogeneity as it would promote species more tolerant or resistant to grazing in communities of intermediate or low preference to the herbivores (Pakeman and Nolan 2009), increasing their similarity to communities of higher preference. Turnover in time (temporal beta diversity) should be higher in intermediate and lower preference communities as grazing increases, whereas turnover in space (spatial beta diversity) should reduce under increased grazing (Salgado-Luarte et al. 2019). Reduced grazing should allow for the increased abundance of grazing sensitive species within preferred and intermediate communities, so turnover in time should be higher in preferred communities but turnover in space should be reduced. If intermediate levels of grazing are continued, as in this study, turnover in both space and time should be lower than under situations where grazing increases or decreases. However, the rate of response of communities to changed grazing levels may differ between increased and decreased grazing; a slower response to reduced grazing is expected (Olofsson 2006). Also, changes in the identity of grazer can also influence vegetation dynamics; more selective grazers such as sheep reduce diversity compared to less selective grazers such as cattle (Tóth et al. 2018).
Species turnover may or may not lead to changes in ecosystem function (Pakeman and Lewis 2018), depending on the level of functional redundancy (the degree to which the loss of a species affects ecosystem functioning depends on the presence of species with similar functional traits, Biggs et al. 2020). However, many of the response traits associated with grazing, such as leaf dry matter content and specific leaf area (Fortunel et al. 2009;Pakeman 2011) are linked to ecosystem function. For instance, leaf dry matter content is a key driver of litter turnover Wilson et al. 1999) and specific leaf area is correlated to productivity (Poorter and De Jong 1999). Therefore, changes in functional beta diversity (the variability in the species functional characters among different communities, Ricotta and Burrascano 2008) should occur if there are changes in abundance weighted species beta diversity as species respond to grazing. We focussed our analysis on vascular plants and on traits known to be responsive to grazing in the experiment (Pakeman and Fielding 2020), namely community weighted mean canopy height and leaf dry matter content declined at higher grazing intensity, specific leaf area increased at higher grazing intensity and start of flowering became earlier. As in many systems these are key 1 3 response traits that are affected by grazing management (Díaz et al. 2007). However, in this type of system these response are also key effect traits as system canopy height impacts nesting density and habitat provision for many invertebrate species (Evans et al. 2015), leaf dry matter content affects livestock foraging choices and nutrient dynamics (Fortunel et al. 2009;Pakeman 2014), specific leaf area affects ecosystem productivity (Poorter and De Jong 1999) and flowering time affects resource timing for pollinators (Robleño et al. 2018).
Using a large-scale, long-term grazing experiment on a mosaic of upland communities in the southern Highlands of Scotland with a long history of grazing (Evans et al. 2005), the following hypotheses were tested: (H1) Increases or decreases in grazing intensity after low-intensity grazing or changes in the identity of the grazer (from sheep only to mixed sheep and cattle) increase temporal beta-diversity of species and functional traits; (H2) Changes in beta diversity are larger in lower preference communities under an increase in grazing and changes are larger in higher preference communities under reduced grazing; (H3) Increases or decreases in grazing away from low-intensity grazing lead to a homogenisation within the mosaic of communities present. Modern methods of assessing beta diversity separate contributions from nestedness, where species loss or gain causes speciespoor sites to resemble a strict subset of species-rich sites (nestedness resultant dissimilarity, Baselga 2010Baselga , 2012) and substitution of one species by another (replacement, Lennon et al. 2001). As plant abundance data was available from the experiment, species beta diversity (d BC ) can be partitioned into dissimilarity derived from balanced variation in abundance between sites (d BC-bal , also known as turnover) and derived from unidirectional abundance gradients (d BC-gra , nestedness, Baselga 2013). Functional beta diversity (or functional dissimilarity Fun-β) can likewise be partitioned into functional turnover Fun-β turn (communities host different functional strategies) and nestedness-resultant functional dissimilarity Fun-β nest (one community hosts a small subset of the diversified functional strategies present in another, Villéger et al. 2013). Therefore, for each of the above hypotheses the question was asked whether changes in beta diversity were driven by balanced variation or abundance gradients.

Experimental details
The grazing experiment is situated on the Glen Finglas estate (56°16′N, 4°24′W), Southern Highlands of Scotland, UK. The vegetation is a mosaic of upland, unimproved plant communities including upland wet heathland (according to the UK National Vegetation Classification, Rodwell 1991Rodwell , 1992 mainly M15 Scirpus cespitosus (now Trichophorum cespitosum)-Erica tetralix wet heath, the mire communities M6 Carex echinata-Sphagnum recurvum/auriculatum (henceforth Carex mire), M23 Juncus effusus/ acutiflorus-Galium palustre and M25 Molinia caerulea-Potentilla erecta (merged as there was little M23 and together called Molinia mire), and the acid grassland communities U4 Festuca ovina-Agrostis capillaris-Galium saxatile (Agrostis-Festuca grassland), U5 Nardus stricta-G. saxatile (Nardus grassland) and U20 Pteridium aquilinum-G. saxatile (bracken). The vegetation types were distributed relatively evenly across the plots within blocks, with some differences between blocks, e.g., Block F had a higher proportion wet heath whilst Blocks A and B had a high proportion of Molinia mire . The climate of the area is cool and moist; interpolated mean temperatures for January and July were 3.04 and 14.12 °C, respectively, and the mean annual rainfall was 2230 mm (1981-2010, Perry and Hollis 2005. Communities were assessed as preferred-Agrostis-Festuca grassland, moderately preferred-Molinia mire, bracken, Carex mire, and unpreferred-Nardus grassland and wet heath-based on observational studies (Hunter 1962) and from the capacity of different communities to sustain grazing animals from stocking rates in experimental studies from the same region (Pakeman 2004(Pakeman , 2014. The experiment was made up of six blocks in three groups of two (Supplementary Material Figure A1). Each block contained four plots of 3.3 ha, roughly 180 m × 180 m, which were fenced in the winter of 2002-03. One of four grazing treatments were allocated randomly to the plots in each block in spring 2003. The treatments were (i) Continuedgrazing at the same low intensity as previous management at three ewes (Ovis aries) per plot (0.9 ewe ha -1 ) to act as the control for the experiment, (ii) High-a tripling of sheep numbers to 2.7 ewe ha −1 , similar to levels seen on farms focussing on profitability, (iii) Mixed-a partial substitution of sheep by cattle (0.6 sheep ha −1 and two cows (Bos taurus) per plot each with a single, suckling calf grazed for 4 weeks in autumn) keeping the overall offtake the same as Continued to assess the effectiveness of cattle at increasing heterogeneity and diversity, and (iv) None-no livestock present to represent agricultural abandonment. Sheep were removed during the winter (December-April) and for short periods for routine farming operations.

Vegetation sampling
Twenty-five fixed points per plot were marked with a wooden peg on the intersections of a square grid at 40 m intervals (Supplementary Material Figure A2). The shape of the grid was adjusted for plots that deviated from square to maintain a similar coverage of the plot. At each fixed point, a vertical pin-frame of five pins, orientated consistently (downhill) from the fixed point, was used to measure the cover of individual plant species and litter (total 125 pins per plot). Values for beta diversity from point quadrat methods such as this will be relatively high compared to areal-based measures (Podani et al. 2013). The vegetation at the field site is dominated by long-lived perennials, many of which form tussocks and an areal-based quadrat method would see little turnover, but contacts by pins with individual leaves will be more dynamic. Sampling was carried out in 2002 (before the fences were erected), 2005,2008,2011,2014 and 2017 in late July/early August, prior to cattle being introduced to the plots.

Trait information
Functional trait data were assembled from two sources; BiolFlor (Klotz et al. 2002) for start of flowering and LEDA (Kleyer et al. 2008) for canopy height, leaf dry matter content and specific leaf area. The datasets are from northern or north-west Europe respectively, so trait values should be representative for the study and trait values were available for all species. All data were checked for consistency of units and outliers.

Statistical analysis
Temporal species beta diversity between species abundances (needle hits) at each sampling point were calculated between 2002 and 2017 using the function beta.pair.abund in the betapart package (Baselga et al. 2018) in R ver. 3.6.1 (R Core Team 2019) with the index family set to Bray-Curtis (Baselga 2013). Changes through time of the three components of species beta diversity, d BC-bal , d BC-gra and their sum d BC , were then assessed with the following two fixed models, Treatment (H1) and Treatment × Community (H2), with Plot nested with Block as the random term using a linear mixed model in lme4 (Bates et al. 2015).
Within plot species beta diversity was calculated separately for 2002 and 2017 by taking the mean of all possible (300) between point measures of beta diversity for each plot. Changes in the mean through time were tested using a linear mixed model with Treatment × Year as the fixed model (H3) and Plot nested within Block as the random model.
Dimensionality of the trait data was reduced from four to two by running a Principal Coordinates Analysis (PCO) on the species x trait matrix using the function dudi.pco in ade4 (Dray and Dufour 2007). Temporal functional beta diversity between 2002 and 2017 for all sampling points was calculated using the function functional.beta.pair from the betapart package (Baselga et al. 2018) with the input matrices being the proportion of each species found at each point and the PCO species coordinates and index family set to Jaccard. The degree of change of functional turnover Fun-β turn , nestedness-resultant functional dissimilarity Fun-β nest and their sum Fun-β through time were tested with the same linear mixed models as temporal species beta diversity above (H1, H2).
The function functional.beta.pair was then used to calculate within plot functional beta diversity in 2002 and 2017, and the mean of all possible between point pairs per plot tested with the same linear mixed model design as within plot species beta diversity (H3).

Species beta diversity
Temporal species beta diversity (d BC ) was similar across all treatments (H1, Fig. 1, Table 1). However, the balanced variation in abundance (d BC-bal ) through time was higher for the ungrazed treatment (None, p = 0.043) than the Continued treatment. In contrast, d BC-gra (the unidirectional abundance gradient) was lower for all the treatments than for the Continued: High p = 0.042, Mixed p = 0.022 and None p = 0.039. Fig. 1 Temporal species beta diversity of the four treatments Cont-Continued (0·9 ewe ha -1 ), High (2.7 ewe ha −1 ), Mixed (0.6 sheep ha −1 and two cows with calves 4 weeks in autumn) and None (no livestock present). Species beta diversity is divided into d BC-bal (balanced variation in abundance), d BC-gra (abundance gradients) and their sum d BC

3
Temporal species beta diversity was lowest in the Molinia mire community (M25) community (H2, Fig. 2). Both the Agrostis-Festuca (U4) and Nardus (U5) grasslands had higher values of d BC-bal , lower values of d BC-gra and higher overall d BC than wet heath, the comparator against which the others were tested (p values between < 0.001 and 0.005, Table 2). Wet heath was chosen as the comparator as it is likely to change slower than the other communities. Analysis of interactions between treatments and vegetation community in the High treatment showed that for the Agrostis-Festuca grassland both d BC-bal (p = 0.033) and d BC (p = 0.031) were reduced compared to the Continued treatment. In contrast, in the None treatment this community showed a reduced d BC-bal (p = 0.018) compared to Continued but an increased d BC-gra (p = 0.019). Similarly, for the Nardus community d BC-bal (p = 0.001) and d BC (p < 0.001) for the High community were reduced compared to the Continued treatment. There were no differences reported for overall treatment changes or community by treatment interactions for the Carex mire (M6), Molinia mire (M25) or bracken (U20) communities.
The unidirectional abundance gradients (d BC-gra ) part of spatial species beta diversity declined across the whole experiment between 2002 and 2017 (H3, Fig. 3, p = 0.042, Table 3). There was also an overall reduction in d BC-bal (p = 0.013) and d BC (p = 0.041) in the High treatment compared to the Continued treatment.

Functional beta diversity
There was no change through time in functional beta diversity (H1, Fun-β, Fig. 4). However, the None treatment saw an increase in functional turnover (Fun-β turn , p = 0.004) and reduced nestedness-resultant functional dissimilarity (Fun-β nest , p = 0.022) compared to the Continued treatment (Table 4).  The absence of an overall change in Fun-β was the same for all communities (H2, no difference from the reference community, wet heath M15, Fig. 5). Analysis of interactions between treatment and community showed only that in the None treatment, the Nardus (U5) community showed a lower Fun-β turn (p = 0.041) and Fun-β (p = 0.005) than the M15 community (Table 5), whilst there was an overall higher level of Fun-β in the Molinia mire (M25) community than in wet heath, but no interaction with treatment.
Spatial functional beta diversity did not change through time across the whole experiment (H3, Fig. 6). However, overall functional beta diversity declined in the High treatment (Fun-β, p = 0.003) driven mostly by a decline in Fun-β turn (p = 0.089, Table 6). There was weaker evidence for this reduction for the Mixed treatment (p = 0.052).

Species beta diversity
Of the two constituent parts of species beta diversity, the balanced variation in abundance greatly exceeded that derived from abundance gradients. However, there were differences between treatments and vegetation types. At the plot level, the higher temporal balanced variation in abundance in the None treatment suggest that species reordering in dominance was highest with the removal of grazing (H1, Fig. 1). In effect, taller growing, grazing sensitive species such as Calluna vulgaris, Deschampsia cespitosa and Vaccinium myrtillus increased whilst there was a decrease in shorter growing species that are tolerant to grazing like Anthoxanthum odoratum and N. stricta . What was surprising was that the d BC-gra (the unidirectional abundance gradient) was lower for all the treatments compared to the continued treatment. The opposite pattern was expected as the  other treatments involved a change in grazing management. However, the experiment has seen increases in height over time, indicating that grazing management is not keeping on top of productivity, and some of the changes in individual species have been highest in the continued treatment, e.g., Narthecium ossifragum and V. myrtillus . The wet heath (M15) and Molinia mire (M25) communities had the lowest turnover (H2, Fig. 2), perhaps because they contain relatively slow growing species, however the wet heath community had the highest unidirectional abundance gradient as it contains species capable of responding to reduced grazing . The Agrostis-Festuca (U4) and Nardus (U5) grasslands had the highest levels of balanced species replacement, perhaps reflecting their more dynamic vegetation compared to the other habitats (Pakeman 2004). Interactions between time and treatment showed that the Agrostis-Festuca (U4) and Nardus (U5) communities had the lowest temporal species beta diversity and, in particular, the lowest replacement of species in the High treatment compared to the Continued. Effectively, only the high treatment, with its increase in grazing, is more closely maintaining the original species composition than the other treatments. The None treatment in the Agrostis-Festuca grassland (U4) community also saw lower species replacement but increased nestedness, reflecting species losses as the taller species increase in cover at the expense of the shorter. This runs counter to the results from principal response curves (PRC) of compositional change for this community , which showed minimal compositional change compared to the Continued treatment. It appears that the analysis of species beta diversity reveals a different dimension compared to compositional changes and that changes are occurring in the Agrostis-Festuca grassland (U4) community as a result of the removal of grazing. Species are being lost despite little movement in ordination space compared to the Continued treatment.
The changes in within plot species beta diversity through time (H3) show that the High treatment had the lowest levels of turnover (balanced species replacement, Fig. 3), adding further evidence that the other treatments are all seeing dynamics driven by reduced grazing levels and increasing vegetation heights (Speed et al. 2013).

Functional beta diversity
Assessing changes in functional diversity at the plot level through time (H1) showed that overall beta diversity was similar across treatments, but the None treatment had higher turnover suggesting a replacement of functional strategies as the vegetation showed increased covers of taller and less grazing tolerant species with high leaf investment costs (high leaf dry matter content, Pakeman and Fielding 2020) and reduced nestedness, reflecting a lower difference in functional space occupied between the two dates than the other treatments (Fig. 4).
Of the plant communities, it appears that the Nardus grassland (U5) was the most tolerant of grazing removal as there was less functional turnover (H2, Fig. 5). However, apart from this all communities show very similar total levels of temporal functional diversity.
Across the plot (spatial) functional beta diversity did not change, but this measure of beta diversity declined in the High treatment suggesting a homogenisation of the functional traits displayed by species (H3, Fig. 6). This is driven by an increase in grazing tolerant species with lower stature, reduced leaf investment costs (lower leaf dry matter content) and stoloniferous growth forms (Pakeman and Fielding 2020) in the more grazing sensitive communities such as wet heath under high grazing.

Grazing and diversity in upland vegetation
Three hypotheses were put forward in the introduction. Firstly, did a change in grazing intensity or species of grazer increase temporal beta-diversity of species and functional traits (H1)? Species beta diversity was highest in the None treatment indicating that removal of grazing had the highest impact on species replacement. All the communities are plagioclimax communities dependent on grazing to maintain their open character (Rodwell 1991(Rodwell , 1992 and the None treatment reflects the removal of management that has sustained their open character for thousands of years (Oosthoek 2013;Stevenson and Thompson 1993). However, much of the evidence suggests that the High treatment was closest in maintaining the pre-2002 species complement (Speed et al. 2013), primarily as productivity in the experiment appears to exceed offtake by a considerable margin in the three least grazed treatments even though the Continued and Mixed treatments kept the same level of grazing as was present pre-2002 ). This may be partly driven by historic high levels of nitrogen deposition (Tipping et al. 2019) and, whilst these are reducing, their effects remain in the system for a considerable time (Payne et al. 2019). It may also be partly driven by increases in temperature across the period of the experiment (Brown et al. 2008;Clark et al. 2010) and also reflect the effect of management decisions to reduce stock numbers across the estate made prior to 2002. However, grazing appears to buffer changes driven by other environmental factors such as climate (Speed et al. 2013).
Secondly, did changes in beta diversity differ between vegetation communities (H2). Specifically, were changes larger in lower preference communities under an increase in grazing and were changes larger in higher preference communities under reduced grazing? The evidence from the study showed that the highest preference community (Pakeman and Fielding 2020), Agrostis-Festuca grassland (U4), and one of the lowest, Nardus grassland (U5), showed similar responses to increased grazing in that they had the lowest levels of temporal species beta diversity. This partially supports the hypothesis (Agrostis-Festuca grassland) and partially does not (Nardus grassland). However, the increased nestedness measure for the Agrostis-Festuca community without grazing does support this hypothesis. Considering this in functional terms, the Nardus grassland, one of the two communities assessed as of low preference, was functionally most tolerant to grazing removal, partially supporting the hypothesis. Overall, there was more evidence supporting this hypothesis than against it, but there was no clear separation of communities based on their grazing preference and dynamics suggesting that all vegetation communities are relatively resistant to change and that only a few changes are detectable.
Thirdly, did changes in grazing lead to a homogenisation within the mosaic of communities present (H3)? There was some evidence that supports this in that within plot levels of both species and functional beta diversity were lower in the High treatment than in the other treatments, agreeing with the findings of Tonn et al. (2019) from mesotrophic grasslands and Hodapp et al. (2018) from a multi-site analysis. The Mixed treatment also had a lower functional beta diversity in 2017, though it was less supported (p = 0.052). This suggests that a combination of increased offtake and increased disturbance has led to some homogenisation as more grazing tolerant species have increased in the less preferred communities in these treatments (Salgado-Luarte et al. 2019;Speed et al. 2013).
The results of these analyses add further information to how grazing affects upland habitats in temperate, nutrient poor ecosystems. In this case higher grazing maintains alpha richness ) and the composition of the vegetation recorded at the start 1 3 of the experiment, but at the expense of some homogenisation. Removal of grazing sees some degree of species loss. These patterns can be added to existing information to assess the biodiversity and other benefits of different management similar temperate, nutrient poor systems; patterns may be very different in systems where unpreferred species spread as grazing increases (Augustine and McNaughton 1998). Management decisions for habitat mosaics are not straightforward as evidence from this experiment shows that individual species, the taxonomic diversity of different groups and ecosystem services like carbon storage and agricultural production all show different responses to changes in grazing (Evans et al. 2015;Smith et al. 2014). This shows there is no "single answer" to managing upland vegetation mosaics and, even for a well-studied system, trade-offs have to be made with imperfect information.