Vegetation sampling and data handling

Vegetation sampling was performed in late June and early July on each of the six sampling occasions. Each plot was gridded into 60 (1 m²) quadrats. In 1972, only 30 randomly selected quadrats from a total of 60 per plot were sampled; however, in subsequent sampling years (1974, 1983, 1996, 2013, and 2016), 60 quadrats were monitored. In each quadrat, the vegetation cover was measured using a 1-m long point quadrat-frame with 10 positions and 2 mm diameter pins with four height strata (0-10 cm, 10-20 cm, 20-30 cm, and >30 cm) (Marrs et al., 1986). Point-quadrat frames were positioned using a permanently marked reference system within each quadrat. In 1972, a total of 10 pins were measured at each quadrat, and 5 pins were measured in the subsequent monitoring. At each pin, the number of pin touches was recorded for each species within each stratum (Marrs et al., 1986). At the same time, litter presence was described by counting the pin-touches in the lowest strata (0-10 cm). The raw data can be found at Rose et al. (2020). Subsequently, the collected point data were treated as follows:

1. At each pin, the number of hits was summed for each of the four strata to provide the total species for each individual pin. These data were then converted to presence/absence data using the “decostand” function in the “vegan” package (Oksanen et al., 2019). Taken together, this process reduced all the data collected at a single pin to either 1 or 0 for each species.

2. At each 1 m² quadrat, we calculated the summed presence for each species and obtained an abundance score of between 0 and 5. As the 1972 data were obtained using 10 pins per quadrat, we selected only the 5 pin positions that were sampled throughout the entire sampling period.

Changes in rare species are not independent of changes occurring in the biotic conditions of the entire plant community (Milligan et al., 2016). Therefore, a series of variables related to vascular plants, bryophytes, and lichens were derived to identify whether they also affected the response of S. hirculus. This included three diversity variables (species richness, and both the Shannon–Weiner and Simpson’s diversity indices), six variables based on taxonomy, and a grazing index (Supplementary Table S1). The taxonomic variables defined were vascular plants, mosses, liverworts, herbs, graminoids, and the graminoid family Poaceae, which is an important group of species that can compete with S. hirculus. To define taxonomic variables and litter abundance weighted indices, the summed total of hits per plot and per taxonomic group and litter was calculated (Alday et al., 2021). The grazing index was derived using Critchley’s (2000) suited species score for grazing, where +1 was associated positively with grazing, 0 represented no discernable response to grazing, and –1 was associated negatively with grazing. Nomenclature was standardized using Stace (2019) for vascular plants and Atherton et al. (2010) for bryophytes.

Data analysis

All statistical analyses were performed within the R Statistical Environment (R Core Team, 2018); the “vegan” package was used for all multivariate analysis (Oksanen et al., 2019) and the “mgcv” package for the Generalized Additive Mixed Modelling (GAMM) (Wood, 2006).

First, the plant community composition and successional trends were described using the absence data for each plot. We used a detrended correspondence analysis (DCA) with the “decorana” function and Hellinger-transformed data. Species that were present on only one occasion were excluded from analysis. The DCA ordination axes were correlated with the three diversity measures (species richness and both Shannon–Weiner and Simpson’s diversity indices) and elapsed time (0, +2, +11, +24, +41, and +44 years) using the “envfit” function with 9,999 permutations. As all diversity variables showed similar significant results (P<0.001), only species richness data are discussed here. The plot distributions were then visualized in ordination space using standard-deviational ellipses (SD-ellipses, 95% confidence limits) with the “ordiellipse”function and the centroids of grazing and sheep excluded plots for each year were plotted to illustrate successional trends. In addition, the area of each ellipse within the ordination space was calculated to measure the size of the overall species pool (i.e., β-diversity). To test whether there was a significant relationship between the species community data and the derived variables (diversity, taxonomic, and grazing indices), we used the “adonis” function in “vegan” with 9,999 permutations stratified by grazing treatment (Supplementary Table S1).

Subsequently, these variables were modelled using GAMMs. We chose GAMMs because they do not make a priori assumptions about the functional relationship between variables (Wood, 2006). GAMMs were fitted using the “gam” function to determine the general shapes of the relationships between (a) the abundance of S. hirculus, (b) species richness, (c) litter, (d) vegetation density (i.e., the summing of all species+litter values obtained), and (e) taxonomic and grazing index variables with respect to elapsed time. Here, the dependent variable was fitted with respect to elapsed time, with elapsed time nested within the grazing treatment to account for temporal correlations (random factors). Four knots were used for the regression splines to avoid overfitting the data but still allowing for unimodal or slightly more complex model fits (Wood, 2006). The results of all analyses are presented in Supplementary Table S2, and all the responses discussed are significant (P<0.01) unless explicitly stated.

Community dynamics over time

The DCA analysis produced eigenvalues of 0.411, 0.179, 0.162, and 0.142 and gradient lengths of 3.28, 2.99, 2.66, and 2.60 for the first four axes. The species plot showed a main gradation along Axis 1 (46% variance, Fig. 1A) from a Festuca ovina and Carex demissa community at the positive end, through a community containing Agrostis stolonifera, Anthoxanthum odoratum, Carex dioica, Carex echinata, Eriophorum angustifolium and Juncus articulatus in the mid-axis towards a community of Carex panicea, Cirsium palustre, J. effusus, and Valeriana dioica at the negative end. Axis 2 reflected a gradation from a grassland community with A. odoratum, Agrostis canina and C. panicea with a large bryophyte presence at the negative end through a community of Carex nigra, C. palustre, Rhizomnium punctatum and V. dioica at the positive end (20% of variance, Fig. 1A). The species richness trajectory was correlated positively with axis 1 and correlated negatively with axis 2, that is, towards the bottom-right quadrant (r²=0.32, P<0.001, Fig. 1A).

The ungrazed plot area occupied was larger than the grazed plot, covering it completely (relative areas of the standard error ellipses: ungrazed=7.66, grazed=2.65). A larger area indicated that there was either larger β-diversity or greater heterogeneity in the ungrazed plot than in the grazed plot (Fig. 1B). The temporal trajectories moved in the same negative direction in both plots along Axis 1 towards C. panicea, C. palustre, J. effusus and V. dioica communities, with the ungrazed plot moving further than the grazed plot (Fig. 1C). Both successional trajectories indicated a reduction in species richness until 2013; however, a downward trend was again noted in the latest 2016 sampling (Fig. 1C), indicating a recent small recovery in species richness, albeit with a different species complement. No Sphagnum spp. were detected.

Changes in abundance of S. hirculus, species richness and litter

A relatively low abundance of S. hirculus was noted throughout the experiment. At the start of the study, the grazed plot contained approximately twice the abundance of the ungrazed plot (Fig. 2A), however the abundance declined to approximately zero in both plots within the first 25 years, and it was below the level of detection of the sampling method used between 1996 and 2013 (Fig. 2A). Thereafter, a slight recovery was observed (Fig. 2A). These temporal responses mirrored the changes in overall species richness (Fig. 2B), which showed a decline over the first 25 years (1972-2000), followed by a recovery, albeit not to the initial values (Fig. 2B). When the ungrazed plot was first fenced, there was little difference in species richness between the grazed and ungrazed plots (with approximately 19 species each); however, there was a markedly greater decline in the ungrazed plot than in the grazed plot in the first 25 years, with a minimum of 7 and 9 species, respectively (Fig. 2B). Both the grazed and ungrazed plots partially recovered and attained a similar species richness at the end of the monitoring period (approximately 13 species by 2016, Fig. 2B).

There was a linear increase in the amount of plant litter in both treatments, with a marginally faster increase at the end of the study period in the ungrazed plot than in the grazed plot (Fig. 2C). The grazing index (Fig. 2D) showed a similarly shaped response to that of species richness (Fig. 2D), although the decline and recovery were less marked. Moreover, no significant differences were observed between the grazed and ungrazed plots.

Changes in litter, taxonomic groups and grazing index

Mosses, liverworts, graminoids, Poaceae and herbs all followed a similar response to that of species richness, with an initial reduction to minimum values in the middle of the study period, followed by a recovery (Fig. 3). However, there was a gradation in the reduction, and vascular plants, mosses, and liverworts were the least affected, intermediate, and most affected, respectively (Fig. 3A-C). There was no difference in the response of liverworts between grazing treatments, and of all the groups, it underwent the greatest reduction and experienced the lowest recovery (Fig. 3C). There were minor differences between the treatment responses over time for the vascular plants and herbs (Fig. 3A, F); however, for graminoids and Poaceae, their reductions in the middle of the study period were more pronounced in the ungrazed plots than in the grazed plots (Fig. 3E). The mosses, liverwort and graminoids/Poaceae (ungrazed only) did not recover to their initial conditions (Fig. 3B, C, E).

S. hirculus changes over time

The main focus of our study was to monitor the presence of the vulnerable species, S. hirculus and assess whether the removal of grazing would help in maintaining their population levels. We initially thought that high grazing pressure was a major factor limiting this species, and it was therefore rather surprising to find that the same temporal response occurred in both treatment plots; that is, an initial decline in the species reaching approximately zero after 25 years followed by a slight recovery. The reduction in the early part of the study suggests that these opposing conservation strategies were neither effective nor harmful (Pullin & Knight, 2009). The minor recovery noted recently is encouraging in the grazed plot, and it suggests that the current grazing level is, at best, not harmful. Determining the extent of the recovery in the ungrazed plot is confounded by the fact that sheep grazing has been allowed since 2013, and it is impossible to determine whether this recovery would have occurred without the re-introduction of grazing. Venterink and Vittoz (2008) studied the last remaining population of S. hirculus in Switzerland and argued that grazing and mowing were needed to remove N from the system and maintain a low above-ground biomass of between 152 g m^-2 and 231 g m^-2 and a N:P ratio in the vegetation of between 7 and 8. In 1983, the respective ranges of the living above-ground biomass at our study site were close to these values, at 181 g m^-2 and 245 g m^-2 in the grazed and ungrazed plots, respectively (Marrs et al., 1989). However, there was a large litter component, and the N:P ratios were also much larger (litter mass, grazed=310±32 versus ungrazed=745±94 g m^-2; N:P ratios, grazed=12.1 versus ungrazed=14.1; Marrs et al., 1989). It is thus likely that during the first 25 years, the above-ground vegetation mass and its N load may have been too high to maintain large numbers of S. hirculus.

It is evident that we do not have enough knowledge about the autecology and population ecology of S. hirculus; that is, we have limited knowledge about the environmental factors that control its recruitment, persistence, and mortality (Hornemann et al., 2012; Pfeifer et al., 2006; Wells et al., 1998; among others). We also have no information about the size and persistence of the seed bank in British upland communities, which limits our reliance on recovery from seed as part of an ecological restoration strategy (Klaus et al., 2018; Rowell et al., 1982). Therefore, several important factors reduce our ability to guide the recovery of S. hirculus. Consequently, for more effective in situ conservation management, it is necessary to conduct new studies on S. hirculus autoecology and its population ecology, and to combine such studies with ecological restoration and conservation programs. Singh et al. (2019) highlighted that non-Sphagnum (“brown mosses”) may play an important role in favoring the persistence and survival of S. hirculus as a component of plant communities in ferns. On our site, no Sphagnum species were detected, but three “brown Mosses” (Calliergonella cuspidata, Campylium stellatum, and Palustriella commutata) were plotted close to the community starting positions in the DCA ordination (Fig. 2), although these have since declined.

Community dynamics through time

It was also clear from our results that despite slight differences in composition between the grazed and ungrazed treatments, the successional processes in both plots followed the same trajectory, but at different rates. It is worth recalling that when the plots were established, the sheep stocking levels had just been reduced from 4.4 ha^-1 to 2.0 ha^-1, so there was likely to be some “adjustment” of species abundance in both grazed plots. It was also expected that there would have been some convergence of the grazed plot towards the ungrazed plot following further stock reductions implemented after 2001, and convergence of the ungrazed plot towards the grazed plot following the re-introduction of grazing to the ungrazed plot in 2013. Irrespective of these events, the community composition trends over the last 40 years have been mainly caused by fluctuations in the abundance of dominant species coupled with reductions in sub-dominant species, including some bryophytes. In both treatments, the dominant species balance changed over time from C. demissa and E. angustifolium with fine-leaved grasses towards Geum rivale (ungrazed only) and coarse-leaved graminoids (including Nardus stricta and J. articulatus) predominantly on the grazed plot, and C. nigra and C. panicea on both the grazed and ungrazed plots. These results are essentially similar to those reported elsewhere for the wider British uplands (Mitchell et al., 2017; Ross et al., 2012). The increased abundance of G. rivale reported here has been noted in a similar ungrazed plot on base-rich soil at Knock Fell, also located at Moor House (Marrs et al., 2020; Rawes, 1981). Evidence from the bryophyte life histories also mirrors the large declines previously reported across the uplands (Britton et al., 2018), although the results of this study show evidence of limited bryophyte recovery.

The compositional changes in vegetation noted here in the grazed and ungrazed plots are in agreement with the changes noted in similar experiments in other communities within the Moor House reserve and in the wider countryside. Therefore, it is likely that compositional changes were mainly driven by landscape drivers, such as the atmospheric deposition of SO₂ and NO_x, which have impinged on species diversity elsewhere in the uplands (Pakeman et al., 2019) and on flora across the entire United Kingdom (Smart et al., 2003a; 2003b; 2006). These results are unsurprising, because it is well known that adjacent sites subjected to similar abiotic pressures and sharing the local species pool produce similar successional processes, which makes large successional differences unlikely (Falster et al., 2017). However, this means that managing the grazing had little effect on plant communities, as similar long-term responses were observed in both the grazed and ungrazed plots, thus showing no positive effect on S. hirculus conservation. In addition, the initial compositional changes in the grazed and ungrazed plots did not favor the conservation of S. hirculus, although the last compositional changes from 2013 showed a slight increase in S. hirculus in both plots.

Change in litter, taxonomic groups, and grazing index

An examination of the litter showed a greater increase in the amount of litter in the last few years in both the grazed and ungrazed plots, although it was generally stable over the first 25 years. These results correspond with the conclusions of Mitchell et al. (2018) from an analysis of repeated vegetation surveys in Scotland, who found that the climate has warmed slightly and become slightly wetter (Holden & Rose, 2011). This could favor the development of denser vegetation and litter cover, and it suggests that decomposition processes have slowed over the period (Marrs et al., 1980; 1989). The future impacts of an increase in litter on S. hirculus are currently unknown, but our results from 2013 showed a parallel increase in litter and S. hirculus.

Overall, the taxonomic group responses were similar to those reported for species richness in both treatments, with an initial reduction and recovery in the previous years. However, there was a gradation in species reduction and liverworts were the worst affected group, followed by mosses and then vascular plants. These results show that community drivers other than grazing have occurred, and these have particularly influenced the persistence of bryophytes and liverworts. Recent studies have suggested that SO₂ and N deposition was a major cause in the reduction of bryophyte and liverwort in the last part of the 20th century (Rose et al., 2016). However, SO₂ acidification peaked in the 1970s, and emissions were reduced by 90% by 2008, while N emissions were also reduced over this period, but only by 50% (Monteith et al., 2016; RoTAP, 2012). These reductions in SO₂ and N deposition have been linked to the start of a reversal in species decline, especially among bryophytes and liverworts (Mitchell et al., 2018; Rose et al., 2016). The results of the present study appear to follow these national trends, in that conditions are becoming warmer and there are reductions in the input of atmospheric pollutants, which is affecting the vegetation composition (Holden & Rose 2011; Monteith et al., 2016). However, we accept that such conclusions are only speculative, and further work is needed to unravel the mechanisms involved. Nevertheless, the decreasing trends in the liverwort and bryophyte populations resemble those of S. hirculus and are in agreement with the results of Singh et al. (2019), who showed a relationship between the persistence and survival of S. hirculus and bryophytes.

Interestingly, there were slight differences between the trends of vascular plants and grasses over time with respect to their responses to treatment, and this suggests that grazing is not as important as expected for vascular plants. In addition, the abundances of graminoids and Poaceae, which are favored by grazers in these areas (Alday et al., 2021), were lower values in the ungrazed plots than in the grazed plots. Surprisingly, there was no difference in the grazing index between the grazed and ungrazed plots, and this could be attributed to the low grazing pressure which was not sufficient to impact species other than the specialized grazing species (Hodgson & Illius, 1998). However, the abundance of mosses, liverworts, and graminoids did not recover to initial levels, even after 44 years, and this trend was similar to that of S. hirculus, which suggests that it takes longer to recover populations or community structures than it takes to degrade them (Alday et al., 2013).

We aimed to describe the results of an early attempt to use an adaptive management approach to maintain populations of S. hirculus at the Moor House NNR in the northern Pennines. We must consider that S. hirculus is a species of conservation concern (BSBI, 2020; JNCC, 2019a; 2019b), which its presence was maintained over time in the flush studied. However, its population declined over the 47 years studied, and it reached a minimum value between the years 1995 and 2010, when it neared extinction, but it has subsequently slightly recovered. This response of S. hirculus appears to have been independent of grazing management approaches (grazing vs. ungrazing plots), and it mirrored wider changes in vegetation within the flush and changes reported elsewhere in the wider uplands that are mainly related to reductions in the deposition of N and sulfur (landscape drivers). The long-term results imply that using a simple adaptive management approach of “just managing” the sheep grazing appears to be ineffectual, and it seems that recent changes in atmospheric deposition or climate warning may have helped the observed recent recovery. However, although the recent recovery of S. hirculus brings hope, these results need to be treated with caution, as the current medium-term responses are unknown.