Quantifying Long-Term Urban Grassland Dynamics: Biotic Homogenization and Extinction Debts

: Sustainable urban nature conservation calls for a rethinking of conventional approaches. Traditionally, conservationists have not incorporated the history of the landscape in management strategies. This study shows that extant vegetation patterns are correlated to past landscapes indicating potential extinction debts. We calculated urban landscape measures for seven time periods (1938–2019) and correlated it to three vegetation sampling events (1995, 2012, 2019) using GLM models. We also tested whether urban vegetation was homogenizing. Our results indicated that urban vegetation in our study area is not currently homogenizing but that indigenous forb species richness is declining significantly. Furthermore, long-term studies are essential as the time lags identified for different vegetation sampling periods changed as well as the drivers best predicting these changes. Understanding these dynamics are critical to ensuring sustainable conservation of urban vegetation for future citizens. species of the multidimensional


Introduction
What drives species composition in cities? Aronson et al. hypothesize that the following hierarchical filters drive novel urban ecosystems: "regional climatic and biogeographical factors; human facilitation; urban form and development history; socioeconomic and cultural factors; and species interactions" [1]. Urbanization has pervasive effects on all parts of the environment and directly affects biotic species composition and richness [2,3]. One such effect is that of biotic homogenization [4,5]. Biotic homogenization is the "process by which regionally distinct, native communities are gradually replaced by range-expanding, cosmopolitan, non-native communities" [4]. Recent evidence for biotic homogenization of urban vegetation includes Melbourne [6] and Montréal [7]. Long-term effects of urban ecosystem filters can lead to either homogenization or differentiation (decreased community similarity) through several mechanisms resulting in various potential outcomes [8]. The differential outcomes highlight the complexity and scale dependence of this process [7,[9][10][11]. Why should we care about biotic homogenization? McKinney emphasizes that it is a major challenge to conservation due mainly to its role in the loss of native species richness in local environments and its effect on human perceptions of nature [5]. Familiarity only with greatly altered urban environments, rich in exotic species, desensitizes urban residents to care for the conservation of local native species [5,12].
Reassuringly, in a recent global analysis of 110 cities, researchers found that the majority of the plant species found in urban areas are still native, but that the densities of these species are far less than those found in rural areas [13]. This analysis makes a strong case for urban nature conservation initiatives [13]. However, biotic homogenization is not the only major concern for native urban biodiversity. Another important often underestimated factor influencing biodiversity is the history

Study Area
The study area is situated in the urban and peri-urban areas of the city of Potchefstroom (26°42'53'' S, 27°05'49'' E, elevation 1350 m) in the North-West Province of South Africa (Figure 1). The city has approximately 250 000 residents and covers a 55 km 2 area. The mean annual rainfall is 550 mm with average temperatures ranging between 0°C in winter to 30°C in summer. This paper reports on the third round of vegetation sampling (in 2019) undertaken at the same sampling locations. The previous sampling events took place in 1995 and 2012. The city was established at its present location in 1841 [43]. The urban development history of Potchefstroom is discussed in du Toit et al. [16]. The spatial expansion of the city is shown in Figure 1. Since 2010, urban expansion noticeably increased in the hills and ridges in the study area, encroaching on many woody grassland communities ( Figure 1, the western part of the study area indicated in black). Potchefstroom is situated in the Grassland Biome at the junction of three vegetation units (Andesite Mountain Bushveld, Carletonville Dolomite Grassland, and the Rand Highveld Grassland) [44].
A recent national list of threatened terrestrial ecosystems lists the Rand Highveld Grassland as vulnerable, meaning that the vegetation has "a high risk of undergoing significant degradation of ecological structure, function or composition as a result of human intervention" [45]. In 2009, 50% of the vegetation was already transformed, with only 2% formally protected [46].

Vegetation Sampling
This study used both vegetation data collected in 1995 and in 2012 [16] as well as new resampled data in 2019. The 1995 surveys were part of an urban plant community study [47]. Vegetation sampling was carried out using phytosociological relevés calculating Braun-Blanquet cover-abundance values for each species present in a sample plot [48]. Sampling was done in open grassland areas (Carletonville Dolomite Grassland and the Rand Highveld Grassland) and woody communities. In the study area, woody communities (dominant taxa, trees and large shrubs) form part of the Andesite Mountain Bushveld vegetation type and are located on the hills and ridges in the region. These hills and ridges form habitat 'islands' in the grassland matrix. Sample plot sizes were 4 x 4 m for open grassland plots and 10 x 10 m for woody grassland plots [49]. Of the 43 plots sampled in 1995 and 2012, 40 could be resampled in 2019 representing 24 woody plots and 16 open grassland plots. Two woody plots were lost. In one plot along the main road, all the trees were removed and the other plot was inaccessible due to illegal dumping of the garden and household refuse and rubble. The only grassland plot lost was due to urban development. In each sample plot, we recorded all the species found within the plot. To allow comparison with the original data we converted the 1995 Braun-Blanquet cover-abundance data to mean percentage values [50] and recorded the percentage cover of each species in the 2012 and 2019 plots. We also recorded any additional species immediately surrounding the plot. We compiled two datasets for this study. The first set contained all the species (annual and perennial) recorded in the plots. The second data set contained only the perennial species recorded [32]. The first dataset was used to test for biotic homogenization in species composition over time [6]. The second perennial dataset was used to determine whether the time lags and drivers have changed since 1995 and 2012. In both datasets, we determined the growth form composition in each plot and calculated indigenous (native) species richness (ISR). In the first dataset, we also calculated exotic species richness, and in the perennial dataset, we additionally calculated the disturbance species richness (DSR) [16]. The DSR represents 31 species (indigenous and exotic) typically found in local grasslands in poor conditions due to overgrazing, trampling, bush encroachment and other disturbances [47,[51][52][53][54]. These species are often pioneer species, weeds and invaders (see supplementary information Table 1 in du Toit et al. [16] for the list of identified species). The DSR was calculated to see whether these species richness patterns react differently to landscape changes than ISR patterns.

Urban Landscape Measures
Calculation of landscape context follows the methodology described in du Toit et al. [16]. The urban landscape measures were quantified for 500 m buffer areas around each of the 40 sampling plots. The six urban landscape measures were age of urbanization (AGE), altitude (ALT), road network density of natural areas (RNDN), percentage natural area (PN), density of dwellings (CD) and landscape diversity (H') [16]. RNDN is a new adapted measure described in du Toit et al. [16], which calculates the density of footpaths and roads inside remnant natural areas. This measure thus quantifies the state of fragmentation as an indication of the quality of the natural area. The description of the other measures and its calculation is also explained in du Toit et al. [16]. The landscape measures were quantified for six time periods in the previous study, 1938, 1961, 1970, 1994, 1999, 2006, and 2010. In this study, all landscape measures of the previous six times were used together with newly calculated measures for 2019 to determine if time lags changed for current vegetation patterns. Up to date Google Earth imagery (2019) was imported into ArcMap 10 [55] to calculate the six measures for the current landscape. As in the previous study, the add-in Hawth's Analysis Tools version 3.27 [56] and Primer 6 [57] were used to calculate RNDN, CD and H'.

Data Analysis
To test for biotic homogenization (all species) we used the Bray-Curtis similarity index [4], instead of the commonly used Jaccard distance similarity matrix as in Zeeman et al. [6], because we had abundance data for all time periods [58]. The similarity index was calculated using the 'similarity percentages' (SIMPER) routine in PRIMER 6. SIMPER calculates the closeness of samples within a group based on the percentage contribution of each species [57]. Data were square-root transformed to allow a greater contribution from the rarer species [59]. The Bray-Curtis similarity matrix was also used to ordinate species composition of the different time periods using non-metric multidimensional scaling (NMDS) in PRIMER 6. We used the same similarity matrix to test for significant differences between species composition of the three time periods by performing an analysis of similarity (ANOSIM) in PRIMER 6. If the generated R statistic value is near 1, there are significant differences between the groups, whilst values near 0 indicate no differences [60]. To determine how widespread and rare species affected compositional similarity over time we calculated the same four categories of indigenous and exotic species (widespread indigenous species, rare indigenous species, widespread exotic species and rare exotic species) as in Zeeman et al. [6]. The Gaston's quartile criterion [61] was used to group species into the four categories. This was done by separately ranking indigenous and exotic species based on occurrence per plot and then using the top 25% of the most abundant species and lowest 25% of the least abundant species in each species list as widespread and rare species respectively [6]. This was done for all three sampling periods. The categorized groupings were then used to calculate the percentage of the total occurrence represented by each category in each time period [6]. The percentage of occurrences were then compared to each other.
To measure community dynamics with time, we calculated species turnover and mean rank shifts between 1995, 2012 and 2019 (see Hallett et al. [62]). We compiled a dataset with the percentage occurrence of each species within each sampling period. Species turnover within each sampling period through time was calculated using the ``codyn'' package in R [62]. This analysis calculates the total turnover as well as indicating the proportion of species appearing and disappearing between two time periods. The mean rank shift indicates the degree of species reordering between time periods. High ranks shift values indicate more reordering of species [63].
To determine whether there were significant differences between the ISR and DSR values of different time periods (perennial species only) we performed an ANOVA in STATISTICA (version 13.5.0.17, https://www.tibco.com/). An ANOVA analysis was also performed on perennial forb and grass species richness of each habitat type for the different time periods.
To test for changes in potential extinction debts we combined time series data with comparisons of species richness patterns to current and past landscape variables [22,26]. Generalized linear modeling (GLM) analysis was performed to calculate the best fit model for the ISR and DSR patterns for the open and woody grassland communities [16]. The ISR and DSR values were the respective response variables, with the predictor variables the six urban landscape measures. The exact procedure followed in R is described in du Toit et al. [16]. The results for the 1995 and 2012 models were used unchanged [16]. We performed new analyses for the 2019 ISR and DSR patterns on the 2019 perennial species dataset.

Little Evidence for Biotic Homogenization in Vegetation Communities from 1995 to 2019
Results of the test for biotic homogenization indicated that no significant changes took place from 1995 to 2019. Results of the Bray-Curtis similarity index comparisons between plots of the same year showed that all plots had low similarity (Table 1). There was a slight trend suggesting that the open grassland plots were decreasing in similarity indicating biotic differentiation. The open grassland plots had very low exotic species similarity, but this similarity was explained by the fact that almost half of the plots did not contain any exotic species. Woody grassland exotic species similarity decreased drastically but could be explained by the slight increase in exotic species richness of these plots with time (from 23 in 1995 to 31 and 30 in 2012 and 2019, respectively).  Figure 2). Moreover, results of the ANOSIM analyses indicated that there were no significant changes in species composition between the respective time periods for any of the groupings of species (i.e., all, indigenous, exotic), for both open and woody grasslands (Table 1). Even though the woody grassland R values were significant at p < 0.05 the values were not near 1 but near 0 which showed that there were no real differences between the species composition of the different time periods. Analysis of the percentage occurrence of the widespread and rare indigenous and exotic species also demonstrated that proportionate distributions of these taxa stayed the same with time ( Figure  3). Both the woody widespread indigenous and exotic species showed a slight increase for 2019. Open grasslands, on the other hand, showed a slight decrease in the proportion of widespread indigenous species. To determine whether the suite of species remained the same, we calculated the amount of shared and unique species found at the different time periods. In the case of the open grasslands, 66% of the species present in 2019 were shared between 2012 and 1995, representing 51% of the original species. Thirty-one percent of the original species were not recorded again in the resampled plots. Patterns for the woody grasslands were similar, with 65% of the species recorded in 2019 shared between all survey periods representing 58% of the original species with 26% unique to 1995.

Vegetation Patterns Changed Slightly between 1995 and 2019.
None of the indigenous species sampled in the study area in the past or current sampling periods were vulnerable or endangered (See Supplementary Table S2 Table S1). None of the changes in perennial species richness between the three sampling years were statistically significant, except for forb species richness in open grassland habitats (Table 2).

Time Lags Changed in Both Open and Woody Grassland Communities
Results of the GLM analyses revealed that time lags have changed since the previous sampling period (Table 3). Open grassland ISR now best correlated with the 2010 urban landscape measures, indicting a 9-year time lag. Open grassland DSR also indicated a 9-year time lag, which was longer than for the 2012 sampling period but shorter than for the 1995 sampling period. Woody grassland plots also had changed time lags where the ISR model correlated best to the 1970 landscape, indicating a 49-year time lag. Woody DSR patterns corresponded to a much shorter time lag period than previous sampling periods, now resembling the open grassland models with a 9-year time lag (Table 3). Hanski and Ovaskainen proposed that another signature for extinction debt is the overabundance of rare species [31]. Distribution of species occupancy histograms for perennial species only (Supplementary Figure S1) showed that the community distributions were unimodal [64] with a distinct peak of highest species richness of species recorded in only one site. This in itself is not a conclusive result for overabundance as much evidence exists to indicate that a high species richness of rare species is not unusual in ecology [64,65]. In the case of our local grasslands, due to a lack of baseline information, we do not know if the rare species were more abundant than it ought to be. In 1995, 27% of the species recorded were only found in one plot, this changed slightly to 26% in 2012 and 24% in 2019. However, 55% of the rare species found in 1995 were not recorded again. Of these, forbs represented 50%. In 2012, 48% of the rare species were not recorded in 2019 of which 60% were forbs. In 2019, 40 species were found in only one site of which 60% were forbs.

Importance of Landscape Factors Driving these Patterns Changed
The importance of the potential drivers in each model also changed for the different grassland models ( Table 3). The 2019 open grassland ISR drivers corresponded to those that were important in 1995, but the relationship with RNDN became negative, i.e., plots with low RNDN had high ISR patterns. The DSR model also corresponded to the same driver that was important in 1995 namely PN. In the woody grassland models, however, the important drivers changed drastically. The ISR patterns were now best explained by the negative relationship with the density of dwellings. In the DSR model, the patterns now corresponded to a positive relationship with both PN and landscape diversity.
Regarding the importance of the landscape measures for the three sampling periods-in explaining ISR patterns, both altitude (positive in open grasslands and negative in woody grasslands) and RNDN (primarily positive) were important in five of the six models with CD important in 4 models. For DSR patterns the most important driver was PN (primarily negative) which was important in four of the six models (Table 3).
Urban landscape measures for each of the eight time periods clearly showed how the urban landscape has changed with time ( Figure 6) (for the full range of measure variation for all the time periods, see Supplementary Table S5)

Discussion
Biotic homogenization has not yet taken place in Potchefstroom, although there is clear evidence for a steady decrease in indigenous forb species over time. Moreover, the results of the GLM models showed that time lags in the response of vegetation patterns towards landscape drivers have changed as well as some of the drivers responsible. This demonstrates the importance of long-term studies and the care that should be taken with interpreting results from limited time frame research. Consequently, it is essential to consider the history of the landscape in conservation efforts.

Little Evidence for Biotic Homogenization in Vegetation Communities from 1995 to 2019
Our results showed that, contrary to results from Melbourne [6], grasslands in the city of Potchefstroom are not homogenizing. A few possibilities might explain this, firstly, that the intensity of urbanization is much less than in the highly urbanized metropolitan Melbourne. In a global sense, Potchefstroom can be seen as a small city. Secondly, Potchefstroom is also less than 200 years old with large areas of grassland left relatively intact until major urban expansion started in 1995 [16]. (Table 1). Percentage occurrence results indicated that the widespread indigenous species were decreasing ( Figure 3). Lastly, the lack of homogenization can be due to the low occurrence of exotic species in the study area. Moreover, the fact that the NMDS and ANOSIM analyses showed no significant differences between years revealed that despite high turnover and rank shifts from 1995-2012 a large suite of widespread species did not change much.

Vegetation Patterns Changed Slightly between 1995 and 2019.
The most important trend in the vegetation patterns is the proportionate composition of forb species, especially perennials, which continually decreased with increases in graminoid species. Moreover, the ANOVA indicated that mean forb species richness was significantly lower in open grassland plots for both 2012 and 2019 compared to 1995. The loss of forb species, especially from 1995 to 2012, might be part of the extinction debt that is being paid. Drivers responsible for this trend might also be linked to altered disturbance regimes of fire and grazing. When herbivores were excluded in natural savannahs, taller fast-growing grasses shaded out lower growing forbs causing reduced herbaceous species richness [66]. Another savannah study showed that herbivory is needed to maintain forb species richness [67]. In Australian urban grasslands, it was also found that graminoid species outcompeted other lower growing species, increasing their susceptibility to local extinction and that reduced fire frequencies would potentially exacerbate the situation [32]. The interplay between fire and grazing have been described as follows "fire affects grazing by altering large-scale foraging patterns, and grazers affect fire by reducing fuel loads and altering fire spread in a landscape" [68].
However, in the current study, the last seven years did not indicate an acceleration in the loss of SR but rather stable numbers with slight differences since 2012, with indigenous graminoid species, even increasing in open grasslands (Supplementary Table S1). Decreasing numbers in indigenous species richness in urban areas is in keeping with evidence from cities worldwide, which showed that even though cities still contain a majority of local indigenous species it is at substantially lower densities compared to non-urban levels [13].
Plot level differences, where some plots lost species and others gained species demonstrate that site scale land use and management history are essential in understanding community dynamics [69]. In the current study, some plots that were degraded in 2012 improved in 2019 due to lower human impacts as represented by the reduction in RNDN values in 2019. However, the presence of time lags and the significant loss of indigenous forb richness indicated that despite slight increases in species richness, the evidence for extinction debts reveals that the situation will not improve but will exacerbate in the future if nothing is done about the situation. In their paper on the global synthesis of plant extinction rates in urban areas, Hahs et al. divided cities into three types based on their landscape transformation history, namely cities highly transformed before 1600 AD, those transformed after 1600 AD but before initial floristic surveys were carried out, and those where initial surveys were done before extensive transformation, i.e., intact native vegetation existed during the surveys [29]. Potchefstroom falls into the third category and is, therefore, likely to carry extinction debts due to the difference in species still present and the predicted species richness based on available habitat. Consequently, development was recent enough for native vegetation to still largely persists in remnant patches. However, without conservation efforts, the authors argue that these populations will not be viable in the long term [29].

Time Lags Changed in Both Open and Woody Grassland Communities.
Our GLM analyses showed that the time lags for open grasslands have changed to a 9 year lag for both the ISR and DSR patterns. In terms of the ISR patterns, this could mean that critical habitat loss due to development and fragmentation had been reached in 2010. This might indicate that habitat quality declined to a critical level due to anthropogenic influences [70], loss of grazing or fire suppression. Moreover, due to the known poor seedbank of grassland species [71], the seedbank might be depleted. The similarity of the DSR time lags to that of the open grassland ISR patterns could also indicate that DSR species responded to the same critical threshold as that of the native grassland species. In our previous paper, we hypothesized that the DSR time lags could act as early warning signals [16], the current results reveal that due to the potential critical threshold reached, this seems no longer to be true. Time lags of the woody grassland ISR were still longer than open grassland species, however, the best fit model now showed that time lags were much longer than initially recorded. A possible explanation for these changes might be found in a recent editorial in response to an intense debate on the effects of fragmentation on biodiversity [72]. The authors highlighted a few key points in the debate that included the importance of the temporal scale. They argued that "more habitat fragmentation might result in greater initial changes in some populations and later changes in others. We cannot assume that the initial rates of biodiversity change observed at different levels of fragmentation will remain over longer time periods" [72]. Another factor to consider in explaining changing time lags is the changes in community equilibrium as habitats change in relation to extinction thresholds [31] and as the extinction debt is paid out [26].
More research on global plant extinctions in urban areas should be undertaken because comprehensive surveys exist of historic and recent floral composition in these urban areas [29]. The changing time lags for the three sampling periods demonstrate that we should be careful in interpreting the results of a single survey. Continuous monitoring is essential to understand the underlying patterns and processes. Our study revealed that the history of the landscape affects current vegetation communities and conservationists need to account for history in contemporary observed patterns in nature.

The Importance of Landscape Factors Driving these Patterns Changed
Our results showed that the importance of some of the landscape measures and the direction of the relationships changed between sampling years. This has important ramifications for selecting and using landscape measures in quantifying time lags as well as incorporating the changing dynamics of the strength of drivers that can change over time in the same landscape. Moreover, this again highlights the importance of long-term studies to monitor changes in the landscape. If we want to aid policymakers and conservationists in effectively managing urban vegetation, decision-makers need to know which drivers are responsible for vegetation composition and patterns. However, we acknowledge that drivers we did not include here could potentially be important in explaining the observed patterns. Two such important drivers are grazing and fire as discussed earlier in explaining the changes in vegetation patterns.

Consequently, How will the Presence of Time Lags and Potential Extinction Debts Influence Conservation Strategies and Planning for Resilient and Sustainable Cities?
In Potchefstroom, biotic homogenization has not yet taken place which means that here remnant natural urban vegetation is still worthy of conservation. However, evidence of a significant decrease in forb species richness combined with the calculated time lags reveal that even in Potchefstroom, with no visible biotic homogenization, keeping the conservation and management status quo is not enough. Pro-active management and restoration are essential to ensure the long-term sustainability of grassland diversity. Better management of extant vegetation and the restoration of native vegetation in the landscape can reduce the extinction debt [25,29]. In grasslands the best restoration results were with planting seedlings in disturbed sites, moreover, exotic weed density was also limited by successfully established native species [73]. Other studies showed that success depends on the re-establishment of swards of the correct indigenous perennial grass species [74]. The removal of dominant grass species allows light availability for intertussock forbs [75,76], however, in their study in Australian grasslands Johnson et al. found that seed limitation of forbs was a constraint to the recovery of degraded grasslands as indigenous forb recruitment mostly relied upon seed addition [76]. This finding is also supported by Zamin et al. [77]. Another management strategy that could be employed in fragmented natural areas is the collection of forb seeds from nearby natural areas and sowing it in degraded urban open spaces including along road verges and railway embankments [78]. Such a strategy could have several benefits in terms of improving the environment, recovering degraded areas with native vegetation, increasing biodiversity, creating "ecological continuity" between the inner city and surrounding natural areas, decreasing maintenance, and also increasing awareness and involvement amongst city residents in terms of the importance of protecting native grasslands [78]. However, not all researchers agree that seeding is the best strategy. In dryland ecosystems, seeding often fails because managers do not account for seed dormancy and germination traits [79]. Whereas, in neotropical grasslands, the translocation of topsoil and transplantation of mature forbs and grasses proved to be the most successful restoration technique [80].
Volis and Deng contend that population demographic assessments should be the first step in any conservation project to determine population viability and appropriate population management [81]. Time lags in the response of species to disturbances can obscure the true conservation status of plants [81]. Our study also revealed that drivers previously associated with vegetation patterns changed. This means that revisitation and long-term studies are essential to understand vegetation dynamics and to continuously ensure that the correct drivers are managed and planned for. Moreover, researchers should be careful to attribute causation too readily to drivers gleaned from limited time frame studies, emphasizing the importance of continuous monitoring. As early as 1990, Magnuson pleaded that "in the absence of the temporal context provided by long-term research, serious misjudgments can occur not only in our attempts to understand and predict a change in the world around us but also in our attempts to manage our environment." [21]. Future measures require the mobilization of efforts well beyond current investment in conservation and need to account for both new and old drivers and possible extinction debts.
Supplementary Materials: The following are available online at www.mdpi.com/2071-1050/12/5/1989/s1, Table  S1: Species richness and percentage contributions (annual and perennial) for all comparable sites between the sampling periods (sites that could not be resampled in 2019 were excluded from previous surveys). Table S2: List of all the indigenous species recorded in the comparison sites for all sampling periods. The list includes the conservation status based on the Red List of South African Plants [82] and the number of plots in which the species occurred for each habitat type and sampling period. The bold text and dark grey fill indicate widespread species and the light grey rare species based on Gaston's Quartile Criterion [61]. Plant names according to the Red List [82] and Germishuizen et al. [83]. Table S3 Listed also is the number of plots in which the species occurred for each habitat type and sampling period. The bold text and dark grey fill indicate widespread species and the light grey rare species based on Gaston's Quartile Criterion [61]. Plant names according to Germishuizen et al. [83]., Table S4: Results of the GLM model for 2019. Model coefficients (±SE) are presented, as well as p values. Model intercept (Int), altitude (ALT), age of urbanization (AGE), percentage natural area (PN), road network density of natural areas (RNDN), density of dwellings (CD), Shannon's diversity index of landscape diversity (H'). Akaike's information criterion (AIC) for the model with the highest AIC value (and year), and the model with the lowest AIC (and year, in bold). GLM model coefficients and SE and p values presented here are from models with the lowest AIC values. L2 and L3 represent the second and third levels of a variable classified as a factor., Table S5: List of the range of variation for each calculated urban landscape measure for all time periods. The listed measures are altitude (ALT), age of urbanization (AGE), percentage natural area (PN), road network density of natural areas (RNDN), density of dwellings (CD) and landscape diversity (H'), Figure S1: Distribution of species occupancy histograms for open and woody grasslands for all sampling periods indicating unimodal distributions.