Effects of Environmental Factors on Plant Productivity in the Mountain Grassland of the Mountain Zebra National Park, Eastern Cape, South Africa

: The relationship between plant productivity, measured according to biomass and species richness, is a fundamental focal point in community ecology, as it provides the basis for understanding plant responses or adaptive strategies. Although studies have been conducted on plant biomass and environmental factors, research concerning mountainous grassland areas is scarce. Therefore, the aim of the present study was to examine the inﬂuence of environmental factors on aboveground plant biomass in the mountainous grassland of the Mountain Zebra National Park, South Africa. Biomass distribution was uneven within the park, owing to certain species having relatively higher biomass values. These differences may be attributed to the chemical and physical properties of the soil, including carbon and nitrogen content, soil pH, and soil texture (sand, silt, and coarse fragments). A disc pasture meter was used to collect biomass data. Multiple regression analysis revealed that most environmental factors did not signiﬁcantly inﬂuence plant biomass. The only environmental factor inﬂuencing plant biomass was soil pH; the inﬂuences of other factors were not statistically signiﬁcant. The results of this study elucidate the interactions of environmental factors with plant biomass. Future research could investigate how environmental factors inﬂuence plant biomass, both below and above the ground in mountainous grassland.


Introduction
Vegetation biomass distribution is a vital concept in plant life history, providing the basis for our understanding of plants' responses or adaptive strategies [1,2].Plant biomass determines the ability of an ecosystem to acquire energy, thus playing a vital role in shaping the community structure and the function of the ecosystem [3][4][5][6].This is because vegetation biomass is a key contributor to soil organic matter, which can influence greenhouse gas emissions in terrestrial ecosystems [7].Thus, vegetation biomass is considered to have a particular function in the global carbon (C) cycle [8].The general prediction is that terrestrial ecosystems globally will be affected by climate change.This is because some of the anthropogenic impacts of climate change include changes in precipitation, atmospheric C, and plant biomass production [9,10].Plant biomass contributes to food webs and the functioning of the overall ecosystem [11].With an increase in atmospheric C concentrations over the past several decades, studying plant biomass is essential for understanding vegetation dynamics, terrestrial ecosystem C stocks, and their responses to environmental changes [12].Plant biomass represents the amount of total energy stored in a given unit of surface area during an observation period.It is influenced by the availability of limited resources such as water, carbon dioxide, and nitrogen (N) [13]; however, the interactive effects of these environmental factors on plant biomass remain unclear [14].Environmental factors play an important role in C accumulation and contribute to C transformation [15].Environmental factors influence C absorption, thereby influencing plant growth and biomass [16,17].
Grassland ecosystems cover approximately 40% of the Earth's land surface, and they account for 1/10 of global C storage [18].Consequently, it is critical to investigate the plant productivity of these ecosystems to understand their vegetation dynamics, terrestrial ecosystem C stocks, and response to environmental changes [18].A shift in environmental conditions would affect vegetation biomass [19].Soil nutrients influence the C absorption of plants, thereby influencing plant growth and causing increases in biomass [16,17].
Several studies have demonstrated that plant biomass influences the structure of a community and the function of ecosystems and is influenced by environmental factors [12][13][14].However, little is known about the interactive effects among these environmental factors on plant biomass [20,21].The responses of plant biomass to different environmental factors are complex, and research on the interactive effects of potential environmental driving factors (such as the climate, soil properties, and topographic properties) on plant biomass distribution in the Mountain Zebra National Park (MZNP) is scarce.Therefore, in the present study, we assessed plant productivity across the MZNP and examined the influence of environmental factors on plant productivity.We focused on aboveground biomass (AGB), since it can be measured without causing any disturbance to the topsoil, unlike belowground biomass (BGB) measurements, which require the entire plant to be uprooted.
The objectives of the present study were to determine whether biomass is evenly distributed across the park, whether environmental factors play a role in plant biomass production, and whether different plant genera account for different biomass values.We hypothesized that a higher availability of N would enhance high plant productivity because N is one of the basic nutrients that limits plant growth and, subsequently, plant productivity/biomass production [22].

Study Area
This study was conducted in the MZNP (Figure 1).The MZNP is approximately 28,500 ha in area and is located within the Eastern Cape Province of South Africa (32 • 22 47 S and 25 • 47 88 E) [22].The park is covered with three main biomes: Karoo Escarpment Grassland (53%), Nama Karoo (37%), and Albany Thicket (10%) (Figure 1) [23].The area is characterized by a cool and arid climate, with an average annual rainfall of 400 mm [22].

Field Data Collection
The dataset was collected from different vegetation units across the MZNP (Table 1; Figure 2).Table S1 of the Supplementary Materials describes each of the vegetation units, and the biomass distribution sampling points are indicated in Figure S1 of the Supplementary Materials [24].A total of 52 plots were sampled, and biomass was collected in a transect measuring 50 m.A disc pasture meter was used to estimate plant productivity rapidly and non-destructively [25]; this was used at each meter along the transect to measure herbaceous biomass [26,27].Plant species data were collected along the same 50 m transect using the step-point method, as described by Mentis [28].The plants were identified to their species level using field identification guides and recorded.The plants were then grouped according to their genus level to determine if the genus had any influence on biomass.The samples were standardized across all plots using the equalizing effort technique, taking into account that a total of 52 plots were assessed.The GPS coordinates were recorded using a GPS device (Garmin GPSMAP 64s) to enable the plots to be used as baseline data in the future.The following soil properties were examined: N, soil pH, clay content, silt, sand, slope, bulk density, coarse fragments, and cation exchange capacity.

Field Data Collection
The dataset was collected from different vegetation units across the MZNP (Table 1 Figure 2).Table S1 of the Supplementary Materials describes each of the vegetation units and the biomass distribution sampling points are indicated in Figure S1 of the Supple mentary Materials [24].A total of 52 plots were sampled, and biomass was collected in a transect measuring 50 m.A disc pasture meter was used to estimate plant productivity rapidly and non-destructively [25]; this was used at each meter along the transect to meas ure herbaceous biomass [26,27].Plant species data were collected along the same 50 m transect using the step-point method, as described by Mentis [28].The plants were iden tified to their species level using field identification guides and recorded.The plants were then grouped according to their genus level to determine if the genus had any influence on biomass.The samples were standardized across all plots using the equalizing effor technique, taking into account that a total of 52 plots were assessed.The GPS coordinates were recorded using a GPS device (Garmin GPSMAP 64s) to enable the plots to be used as baseline data in the future.The following soil properties were examined: N, soil pH clay content, silt, sand, slope, bulk density, coarse fragments, and cation exchange capac ity.

Statistical Analyses
Species diversity was measured using the Shannon diversity index (SDI) [29].Other measures such as species richness and Simpson indices [30] were also considered, and both showed strong correlations (87.80% and 95.20%, respectively) with the SDI.The observed degree of correlation between the index pairs was expected, as all the indices are variants of the Hill index [31,32].Therefore, the choice of index for determining the species diversity was established without bias.
Given the proximity of the data points to one another, we first checked for autocorrelation using the Moran index test [33], as implemented in the ape package of the R statistical software version 4.3.0[34,35].The test was conducted using the collected biomass data.The distances between the longitude and latitude measurements were converted to weights using the population density adjustment technique [36,37].A p-value of approximately 0.0094 was obtained.This indicates that a statistically significant autocorrelation between the observations exists.Therefore, all tests of association with biomass presented here were solved as generalized least squares regression problems [38] with a spherical correlation structure, as recommended by Beale et al. [39].The correlation structure was

Statistical Analyses
Species diversity was measured using the Shannon diversity index (SDI) [29].Other measures such as species richness and Simpson indices [30] were also considered, and both showed strong correlations (87.80% and 95.20%, respectively) with the SDI.The observed degree of correlation between the index pairs was expected, as all the indices are variants of the Hill index [31,32].Therefore, the choice of index for determining the species diversity was established without bias.
Given the proximity of the data points to one another, we first checked for autocorrelation using the Moran index test [33], as implemented in the ape package of the R statistical software version 4.3.0[34,35].The test was conducted using the collected biomass data.The distances between the longitude and latitude measurements were converted to weights using the population density adjustment technique [36,37].A p-value of approximately 0.0094 was obtained.This indicates that a statistically significant autocorrelation between the observations exists.Therefore, all tests of association with biomass presented here were solved as generalized least squares regression problems [38] with a spherical correlation structure, as recommended by Beale et al. [39].The correlation structure was formed using the latitude and longitude coordinates.Therefore, despite playing such important roles, these coordinates never explicitly appeared in the outputs of any of the analyses.
The generalized least squares (GLS) method was implemented to optimally extract evidence regarding biomass patterns at the MZNP.Two things about how we implemented this method are noteworthy.First, the all-subsets variable selection technique was used to identify the "best" linear combination of the environmental variables (including the SDI) that explained the observed variation in average biomass per plot [40,41].As a result, 32,767 GLS subset models were fitted.The Akaike Information Criterion [42] was used as the selection criterion.Second, an approximate value was obtained for the regression correlation coefficient (R 2 ).This was necessary because GLS models do not return R 2 values.The reported approximation was obtained as the square of the correlation between the observed (y) and the model-predicted ( ŷ) biomass values.Some adjustments were made to the data to minimize bias.Nine VUs were observed in total.However, only one plot was recorded for each of the two units (VU3 and VU11).Therefore, excluding VU1, neighboring VUs were merged to ensure that each class spanned at least 5% of the plots.The resulting VU classes used in subsequent discussions were VU1, VU2-3, VU4-5, VU7-8, and VU9-11.Furthermore, 13 plots share five latitude or longitude values.Therefore, to guarantee that the autocorrelation model was implemented successfully, duplicated values had to be eliminated.As such, random uniform (0.0001, 0.0003) values were added to the duplicated coordinates.To control the false discovery rate resulting from multiple tests, a Benjamini and Hochberg [43] correction was used to adjust all subsequent p-values.Hence, all p-values reported represent the corrected estimates.Of note, in accordance with the convention, a p-value threshold of 5% was adopted for determining the statistical significance of the evidence present in the analyzed data with regard to the hypotheses considered.All the analyses were performed using R software 4.3.2[35].

Results
First, the biomass pattern across all 2594 (50 × 50 + 47 × 2 plots) locations was examined.A graphical illustration of this is shown in Figure 2, and a quantitative summary is presented in Table 1.The distribution was positively skewed, with only a few of the values greater than 15 g/m 2 .Although the maximum biomass recorded was 36 g/m 2 , 75% of the records did not exceed 10 g/m 2 .Therefore, it is not surprising that the data contained significant evidence (p-value of approximately 0.0000) against the hypothesis that biomass was normally distributed across the MZNP.The p-values provided in red indicate the existence of statistically significant evidence that goes against the hypothesis, i.e., that shows the corresponding data to be normally distributed.
Second, the distribution of average biomass as a function of the observed genus was investigated.See Figures 2 and 3 for a visual illustration and a quantitative summary of this, respectively.The distribution was approximately normally distributed, given the concept of the central limit theorem [40].There was insufficient evidence in the observed data (p-value = ~0.4011) to reject the claim that mean biomass, as a function of genera, was normally distributed.It can, however, be inferred that the distribution showed a mild left tail caused by the low average biomass of the Felicia genus.
The average biomass distribution pattern was further investigated in terms of plots and with respect to species (Figures 2 and 3).In both instances, the inferences were similar to those obtained, with respect to the approximate normal distribution of genera and high pvalues that indicated insignificant evidence against the normal distribution claim (Table 2).As seen for the genera, the other pairs of mean biomass distributions had tails.The tail with respect to the plots was mild and caused by the relatively high average biomasses of Plot 3 and Plot 35.However, the tail in the case of the species was heavy on both ends of the distribution.This was because of the particularly high mean biomass recorded for Asparagus striatus and Diospyros austro-africana, as well as the relatively low average biomass recorded for Felicia filifolia.The two species with large mean biomasses had to be temporarily excluded from the data so that the normal distribution hypothesis would not be rejected.The average biomass computed plot-wise contained the least evidence against a normal distribution hypothesis.Therefore, this justified the utilization of the GLS method to subsequently analyze the relationship between average biomass values and different environmental factors (Table 3).

EER REVIEW 6
Figure 3. Mean biomass distribution with respect to the observed genera, plots, and species.The illustration on the left shows the densities (that is, unitless relative frequencies such that the area under each curve sums to one).The right-hand image highlights the magnitude of the average biomass with respect to each of the data points-Genus: decreases from left to right; Plot: increases from bottom to top, in the center; Species: increases from left to right.
The average biomass distribution pattern was further investigated in terms of plots and with respect to species (Figures 2 and 3).In both instances, the inferences were similar to those obtained, with respect to the approximate normal distribution of genera and high p-values that indicated insignificant evidence against the normal distribution claim (Table 2).As seen for the genera, the other pairs of mean biomass distributions had tails.The tail with respect to the plots was mild and caused by the relatively high average biomasses of Plot 3 and Plot 35.However, the tail in the case of the species was heavy on both ends of the distribution.This was because of the particularly high mean biomass recorded for Asparagus striatus and Diospyros austro-africana, as well as the relatively low average biomass recorded for Felicia filifolia.The two species with large mean biomasses had to be temporarily excluded from the data so that the normal distribution hypothesis would not be rejected.The average biomass computed plot-wise contained the least evidence against a normal distribution hypothesis.Therefore, this justified the utilization of the GLS method to subsequently analyze the relationship between average biomass values and different environmental factors (Table 3).Genus / Plot / Species Species Genus Plot Figure 3. Mean biomass distribution with respect to the observed genera, plots, and species.The illustration on the left shows the densities (that is, unitless relative frequencies such that the area under each curve sums to one).The right-hand image highlights the magnitude of the average biomass with respect to each of the data points-Genus: decreases from left to right; Plot: increases from bottom to top, in the center; Species: increases from left to right.Vegetation unit and water pH ("p-value" approximately 0.0385) were the two variables that were each identified as influencing the average biomass (Table 3).An increase in water pH of 0.10 tended to be associated with a decrease in average plot biomass of approximately 0.069 g/m 2 .The average plot biomass seemed to be similar across Vus, except for VU2-3 ("p-value" ≈ 0.0333).Compared with that at VU1, the average biomass tended to be approximately 3.1 g/m 2 higher at VU2-3.The relationship between average biomass and each of the average yearly rainfall and soil organic C values merits further investigation as the significance of these associations only became weakened after the p-value adjustment.Similar post-adjustment evidence weakening was observed for VU9-11.Therefore, despite the identification of a pair of variables that were significantly associated with average biomass, more extensive studying is required with data collected from more plots.
We further sought to understand the optimal linear combination of the observed environmental factors, including the standard deviation index, which explained the recorded variation in average biomass (Table 4).A three-variable model that included aspect, VU, and water pH was identified as the optimal model.It was estimated that the model explained approximately 43.65% of the observed variation in average biomass.Given that VU and water pH have been previously identified to have a significant pairwise relationship with average biomass, their inclusion in the parsimonious model was not surprising.Based on a similar analogy, the presence of aspect in the model was unexpected.A possible explanation for this is that the aspect variable had a confounding effect.This hypothesis was investigated further.Figures 4 and 5 were used to investigate whether the importance of the aspect variable in the optimal biomass model was because of its potential role as a confounder.For the purpose of representation, the aspect variable used to generate the graphs in the figures was adjusted.North-east and north-west were re-classified as north.Similarly, south-east and south-west were re-classified as south.
Ecologies 2023, 4, FOR PEER REVIEW 8 VU and water pH have been previously identified to have a significant pairwise relation ship with average biomass, their inclusion in the parsimonious model was not surprising Based on a similar analogy, the presence of aspect in the model was unexpected.A possi ble explanation for this is that the aspect variable had a confounding effect.This hypoth esis was investigated further.and 5 were used to investigate whether the importance of the aspect varia ble in the optimal biomass model was because of its potential role as a confounder.For the purpose of representation, the aspect variable used to generate the graphs in the fig ures was adjusted.North-east and north-west were re-classified as north.Similarly, south east and south-west were re-classified as south.Figure 4 shows that VU9-11 had the highest average in the eastern and southern aspects, whereas this did not appear in the west.Moreover, VU2-3 was among those with the highest mean across all aspects, except in the east, where it did not appear.Therefore, the aspect played a confounding role in the association between mean biomass and VU.
More stimulating features were evident in Figure 5, where the confounding effects of the aspect on the association between water pH and average biomass were explored.Except for the south, all the aspect categories had extreme water pH values that skewed the distributions leftward.Notably, the extreme values were influential.The directions of the slopes of the regression lines changed from negative to positive for the eastern and northern aspects.For the western aspect, excluding only the extreme values caused the regression slope to increase.Therefore, it is plausible to claim that there was a negative association between water pH and average biomass for the western aspect.Conversely, for the southern aspect, the association between water pH and average biomass appeared to be strictly positive.This contrast validates the suspected confounding attribute of the aspect variable on the pairwise relationship between the water pH and the average biomass.

Discussion
The MZNP is the interface between three biomes, namely, grassland, Nama Karoo, and Albany thicket.The biomass across the park was not normally distributed, owing to the high biomass of the genera Asparagus and Diospyros as well as the relatively low average biomass recorded for F. filifolia.Moreover, D. austro-africana and Asparagus are dominant within the Nama Karoo biome, while F. filifolia dominates the grassland biome within the park.The high-biomass species were visibly dominant in the study area while F. filifolia was not dominant, which may explain their higher biomass contribution.Soil pH was the sole factor investigated in this study that markedly affected plant biomass production.Other soil factors that were examined did not demonstrate significant effects on plant bi- Figure 4 shows that VU9-11 had the highest average in the eastern and southern aspects, whereas this did not appear in the west.Moreover, VU2-3 was among those with the highest mean across all aspects, except in the east, where it did not appear.Therefore, the aspect played a confounding role in the association between mean biomass and VU.
More stimulating features were evident in Figure 5, where the confounding effects of the aspect on the association between water pH and average biomass were explored.Except for the south, all the aspect categories had extreme water pH values that skewed the distributions leftward.Notably, the extreme values were influential.The directions of the slopes of the regression lines changed from negative to positive for the eastern and northern aspects.For the western aspect, excluding only the extreme values caused the regression slope to increase.Therefore, it is plausible to claim that there was a negative association between water pH and average biomass for the western aspect.Conversely, for the southern aspect, the association between water pH and average biomass appeared to be strictly positive.This contrast validates the suspected confounding attribute of the aspect variable on the pairwise relationship between the water pH and the average biomass.

Discussion
The MZNP is the interface between three biomes, namely, grassland, Nama Karoo, and Albany thicket.The biomass across the park was not normally distributed, owing to the high biomass of the genera Asparagus and Diospyros as well as the relatively low average biomass recorded for F. filifolia.Moreover, D. austro-africana and Asparagus are dominant within the Nama Karoo biome, while F. filifolia dominates the grassland biome within the park.The high-biomass species were visibly dominant in the study area while F. filifolia was not dominant, which may explain their higher biomass contribution.Soil pH was the sole factor investigated in this study that markedly affected plant biomass production.Other soil factors that were examined did not demonstrate significant effects on plant biomass production.This is inconsistent with the findings of the studies by Baraloto et al., Sun et al., and Yang et al. [44][45][46], which report significant relationships between biomass and soil properties in grasslands.Furthermore, Sun et al. [45] reported a positive relationship between soil C density and AGB.
Soil N content did not influence plant biomass production in the present study.This is in line with the findings of several studies stating that soil N content has little to no influence on plant biomass production in most terrestrial ecosystems [47,48].Bhandari and Zhang [49] observed a negative relationship between soil N content and plant biomass production, in which a relatively high N content resulted in a relatively low biomass production, while a relatively low N content translated to a relatively high biomass production.A study conducted in an alpine meadow also revealed that N content is a major factor affecting AGB, influencing the species richness [50].
Vegetation biomass is also an important factor [51] as biomass is a major contributor to soil organic matter, which can influence greenhouse gas emissions in the terrestrial ecosystem.Biomass also plays a major role in energy production, which is one of the greatest essential quantitative characteristics for understanding community structure and ecosystem function.
The positive association observed between plant biomass and soil pH in the present study is consistent with the findings of a study conducted by Bhandari and Zhang [49], which states that soil pH and plant biomass production have a positive relationship at the southern aspect [52].Although most of the environmental factors did not influence biomass distribution in the present study, studies that examine environmental factors and plant biomass remain invaluable for an overall understanding of plant biomass distribution in an area.
Regionally, AGB is primarily determined according to altitude, climate, and soil fertility [53,54].However, in the present study, there was no significant relationship between AGB and mean annual precipitation.
The AGB and environmental factors across 52 sites were analyzed.Based on the results, only one of the soil properties, soil pH, influenced the plant biomass.Although soil nutrients are highly important in plant biomass production, we found their influences on AGB and plant to be insignificant.The results of the present study show that various factors explain the spatial distribution of plant biomass abundance, which may also be influenced by various vegetation types, geology, and rainfall of an area.Furthermore, the spatial distribution may be affected differently by different environmental factors.The main limitation of our study was that only AGB was considered, rather than the whole plant biomass, which includes both the above-and belowground biomass.Future studies should focus on the whole plant, as this may provide a better understanding of how soil properties influence plant biomass production and allocation.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/ecologies4040049/s1. Table S1.Management Units at Mountain Zebra National Park and their description.Funding: This work was funded by the South African National Park and the JRS6 Biodiversity Foundation SPKH1400 in Washington, DC, through South African National Parks.The sponsors were involved in the running costs of data collection.The sponsors were not involved in the study design, the data collection and analysis, the writing of the report, or the decision to submit the article for publication.

Ecologies 2023, 4 ,Figure 1 .
Figure 1.Location of the Mountain Zebra National Park study area showing its three biomes and the location of the park within South Africa.

Figure 1 .
Figure 1.Location of the Mountain Zebra National Park study area showing its three biomes and the location of the park within South Africa.

Figure 2 .
Figure 2. Distribution of plant biomass across the 2594 surveyed land plots.The vertical axis shows the densities (unitless) of the observations that are within the respective intervals.In other words, the y-axis is the relative frequency, scaled such that the area under the histogram sums to one.

Figure 2 .
Figure 2. Distribution of plant biomass across the 2594 surveyed land plots.The vertical axis shows the densities (unitless) of the observations that are within the respective intervals.In other words, the y-axis is the relative frequency, scaled such that the area under the histogram sums to one.

Figure 4 .
Figure 4. Conditional distribution of average plot-wise biomass across vegetation units.Bar heights and half-arrow widths are the means and standard errors of the average biomass values, respec tively.

Figure 4 .
Figure 4. Conditional distribution of average plot-wise biomass across vegetation units.Bar heights and half-arrow widths are the means and standard errors of the average biomass values, respectively.

Figure 5 .
Figure5.Exploration of the pairwise relationship between average biomass and water pH, conditional according to the aspect.The plot on the right exhibits the generalized least squares (GLS) regression lines restricted to water pH values between pH 7 and 8 to avoid extrapolation.The dashed lines denote cases where extreme water pH values (<pH 4) were excluded.The unbroken lines were generated from the complete data.

2 )Figure 5 .
Figure5.Exploration of the pairwise relationship between average biomass and water pH, conditional according to the aspect.The plot on the right exhibits the generalized least squares (GLS) regression lines restricted to water pH values between pH 7 and 8 to avoid extrapolation.The dashed lines denote cases where extreme water pH values (<pH 4) were excluded.The unbroken lines were generated from the complete data.
Figure S1.Location of the study area showing biomass distribution across the vegetation units in the Mountain Zebra National Park, South Africa.Author Contributions: Conceptualization, N.M., A.R., S.A. and H.B.; methodology, N.M.; formal analysis, N.M. and H.S.; data curation, N.M.; writing-original draft preparation, N.M.; writing-review and editing, A.R., S.A., H.B. and H.S. All authors have read and agreed to the published version of the manuscript.

Table 1 .
Details of the environmental variables used in the present study.

Table 1 .
Details of the environmental variables used in the present study.

Table 2 .
Biomass distribution as a function of all considered locations, observed genus, species present, and plots.The p-values provided in red indicate the existence of statistically significant evidence against the hypothesis that the corresponding data are normally distributed.

Table 2 .
Biomass distribution as a function of all considered locations, observed genus, species present, and plots.The p-values provided in red indicate the existence of statistically significant evidence against the hypothesis that the corresponding data are normally distributed.

Table 3 .
Independent generalized least square regressions fitted with average biomass per plot as the response.
* The p-values in red highlight statistically significant associations.The underlined p-values indicate the associations that were significant prior to p-value adjustment.

Table 4 .
Analysis of variance for the "best" subset of the generalized least squares (GLS) model.

Table 4 .
Analysis of variance for the "best" subset of the generalized least squares (GLS) model.