Contrasting the Effect of Forest Landscape Condition to the Resilience of Species Diversity in a Human Modified Landscape: Implications for the Conservation of Tree Species

: Using landscape moderation insurance and Intermediate Disturbance Hypothesis (IDH) as frameworks, this study assessed the response of local assemblage among different land use regimes (mean β -diversity), using the Jaccard dissimilarity matrix in contrasting Human Modified Forest Landscapes (HMFLs). The study was conducted at the relatively simplified Mafhela Forest Reserve and the complex Thathe Vondo Forest Reserve in South Africa. The patterns of overall β diversity between HMFL and State-protected Indigenous Forests (SIF) were compared and the leading change drivers were then untangled. This study found that human disturbance affects mean β -diversity of local assemblages among land use regimes between the two HMFLs in an ecologically contrasting manner. The HMFL in Mafhela Forest Reserve had distinct local assemblages among land use regimes and did not conform to the expectation of IDH. On average, HMFL had the same average local species richness as SIF, mainly due to change in species composition (species replacement) induced by land use disturbance. Land use intensity gradient was the leading change driver to explain the overall β -diversity of the Mafhela Forest Reserve. The findings in the Thathe Vondo Forest Reserve were in contrast with the Mafhela Forest Reserve. Although on average the HMFL had the same local species richness as SIFs, this was mainly due to a trade-off of species gain in trees along the rivers and streams and species loss in Culturally Protected Areas (sacred forests) (CPA) as expected by IDH. The contrasting findings imply that the effectiveness of any alternative conservation strategy is context-dependent. The resilience of local assemblages and conservation value of HMFL depends on the condition of the overall forest landscape complexity and cannnot be captured by one theory, nor by one species diversity matrix (e.g., β -diversity or Richness). It thus demands the application of complementary theoretical frameworks and multilevel modeling.


Introduction
It has been argued that protected areas can neither serve as a standalone strategy to protect rare or endangered species [1] nor be effective enough to tackle the global biodiversity loss from anthropogenic disturbance [2]. Also, despite increasing hope for the conservation potential of forests and trees in modified rural landscapes, there has been growing controversy focusing on the richness, composition, and survival of biodiversity given persistent anthropogenic disturbances [3]. Such Secondly, every landscape is a unique socio-ecological system [26]. The spatial patterns of extent, frequency, and intensity of resource extraction in a landscape are driven by elevation [27], accessibility, distance from the villages, and availability of preferred species for specific uses [28]. This suggests that a particular forest landscape complexity emerges out of a myriad of interactions of anthropogenic and environmental change drivers [26] that are superimposed on the original forest biota. Whether the original biota of a landscape remains resilient depends on the spatial patterns of land use intensity, cross-scale interaction of disturbance regimes [29] and the conditions of existing forest complexity [30]. The condition of existing forest complexity includes the condition of forest cover, ecological connectivity (composition and configuration of a species pool) and successional stages of different patches [15].
According to the landscape moderated insurance hypothesis [30], local assemblages in a complex landscape are expected to have better resilience and stability of ecological processes, even under a continuously changing environment. Better ecological connectivity of a complex landscape enhances the supply of propagules from the species pool and stronger neighborhood effects among land use regimes of different intensity during disturbance-recovery dynamics. Hence, natural factors (e.g., elevation, slope gradient, position) obtain major ecological importance when the human influence on pre-existing forest complexity is minor or when human influence is widespread and fairly uniformly distributed across the whole landscape. The converse is true when land use disturbance overrides and simplifies the condition of the pre-existing environmental conditions [31]. In that context, the existing landscape complexity moderates the effect of land use on β-diversity [32] at a different spatial scale [30]. Consequently, many studies have suggested that the potential of the human modified landscape for conservation is context-dependent [3]. However, to the best of our knowledge, there has not been any study that reveals the conditions under which both landscape moderation insurance and Intermediate Disturbance Hypothesis (IDH) provide a reliable predication.
The primary objective of this study was to investigate the effect of land use regimes and environmental variables in the HMFLs on the resilience of species diversity, under different conditions of forest landscape complexity. The study was further intended to reveal the consequence of land use regimes of different intensity gradients on the cumulative contribution of the human modified landscape to the conservation of tree species, in comparison with relatively well-protected forest areas. This was crucial as the few remaining tropical forest fragments are found embedded within the same human dominated biosphere [10].
Drawing upon the above arguments and scientific theories, the study hypothesized (H) the following: H1. The condition of forest landscape complexity determines the pattern of mean β-diversity and the conformity of species richness response to IDH along the land use intensity gradient of human modified forest landscapes (HMFL);

H2. The condition of forest landscape complexity determines the local contribution of the HMFL to species richness and overall β-diversity (LCBD) in comparison with the strictly protected area;
H3. The influence of land use gradient and other environmental change drivers to overall β-diversity of a particular Forest Reserve depends on the conditions of forest landscape complexity.

Study Area
This study was conducted at Mafhela Forest and Thathe Vondo Reserves in the Vhembe Human Biosphere Reserve (VBR) located in Limpopo Province of South Africa. The two Forest Reserve areas (FR) belong to the eastern part of Soutpansberg Mountain Forest complex that stretches from Louis Trichardt to Thohoyandou. Thathe Vondo (TVFR) and Mafhela Forest Reserves (MFR) are located at 22°52′ S, 30°20′ E and 23°01′ S, 30°30.35′ E, respectively ( Figure 1). Both Reserves have an altitudinal range of 700-1700 m above sea level. The areas receive rainfall from October to March (on average-724 mm), and the average temperature ranges from 35 °C in summer to 18 °C in winter. Both Forest Reserves retain some of the remaining moist forest cover in South Africa [33]. Notwithstanding, local people residing within the boundaries of the Reserves still practice various tree-based traditional land use, whose potential impact on conservation is unknown. For this study, the existing land use regimes were classified, in consultation with traditional leaders, based on the configuration of trees, their typical land use, cultural practices, protection gradient, and their management system [34,35]. The existing tree-based traditional land use regimes of the two Forest Reserves were grouped into two major groups as follows: (a) Human modified landscape (HMFL): consisting of three tree-based traditional land use regimes under the custody of traditional authorities and the local community. These are: (b) Trees Along Streams and Rivers (TATR): local community members are not allowed to harvest live trees, but occasionally access the place for livestock grazing, watering, and shading (intermediately disturbed); (c) Common Resource Use Zones (CRUZ): this is an open access area for the harvesting of wild food, construction materials, livestock browsing and grazing, traditional medicines and others (highly disturbed); (d) Culturally Protected Forest Areas (CPA): these include sacred/holy forests that are protected by royal families for cultural values and only accessible to them (minimally disturbed). (e) State indigenous forests (SIF): these are fragmented forest patches, with minimal to no human disturbance and are legally protected by government conservation agencies.

Sampling Design
Three different asymmetrical nested sampling designs with hierarchical factors relevant to the hypotheses and objectives of the research were adopted after Anderson et al. [36]. To analyze the first hypothesis on β-diversity along the land use gradient in HMFL of the FR and the conformity to IDH, two factors -Transect (Tr) and land use regime (La) were considered. The transect was nested under land use regimes. To determine the effect of land use on β-diversity, land use regimes were kept as a fixed variable (Figure 2). To compare the overall species richness and the local contribution to β-diversity (LCBD) of HMFLs against SIF (HMFL Vs. SIF) for the second hypothesis, three factors were considered; Transect, land use, and HMFL Vs. SIF. The transects were nested in land use regimes. Land use regimes were nested in HMFL vs. SIF, and HMFL vs. SIF was kept as a fixed factor ( Figure 3). Lastly, for the third hypothesis to explore the influence of land use regimes and other environmental change drivers to overall β-diversity of a particular forest Reserve, the plot was nested in Forest Reserves (FR). Each FR consisted of three land use regimes of the HMFLs (CRUZ, TATR, and CPA) and SIF. Identification and location of land use regimes within the Forest Reserves were done with the guidance of local informants [37,38]. In each land use regime, five transects (Tr) were established. Each transect was 50 m long and separated from each other at least by 200 m. In each transect, three (3) 20 m × 10 m rectangular plots (P) were established and spaced 10 m apart along a linear transect. However, only three transects were used for data collection in TATR in Mafhela Forest Reserve as a large part of the remaining blocks of forest patches had been cleared for horticultural production. Similarly, only three transects were used for data collection in CPA in Thathe Vondo Forest Reserve. This was due to a restriction of access to the other part of the land use regime by the local chief related with traditional beliefs.
Previous studies in many tropical countries showed that a 40 m × 5 m transect is appropriate for tree species diversity survey [39]. However, in this study, the transect length was 50 m long and each transect was separated from the other ones by at least 200 m. However, unlike other similar studies (e.g., Gillison [40]) that placed four (4) 5 m × 5 m square shaped sample plots in each transect, this study established three (3) 20 m × 10 m rectangular plots (P) that were spaced 10 m apart along a linear transect.
The modification of Gillison [39] and Gillison et al. [40] sampling approaches was important to strike a balance between the observer's efficiency of data collection and effectiveness of the sampling effort. Experiences in different countries showed that observer fatigue increases if a transect size larger than 40 m × 5 m is used in complex vegetation [41]. Gillison [39] discussed the rough train as a likely factor for observer's fatigue, which also characterizes the ragged mountain forests of this study area. Hence, the slight increase in transect length combined with the reduction in the number of sample plots per transect was intended to alleviate observer fatigue. At the same time, rectangular sample plots have been proven to capture more species than square shaped sample plots of the same size [41], let alone when using the relatively larger sample plot size covered in this study. The overall sampling area coverage per transect combined with the slight increase in transect length was an attempt to capture a board range of local habitat heterogeneity.

Measurement of Tree Species Assemblage
In this study, all perennial woody plants with a diameter at breast height (dbh) ≥ 2 cm and height of ≥2 m were considered as trees and enumerated. The scientific and vernacular names (from local informants) of observed tree species in each plot were recorded. In instances where tree species identification was not possible in the field, tree voucher specimens were collected and later identified at the Thohoyandou Botanical Garden and Herbarium.

Change Drivers of Overall β-Diversity of FRs
Land use regime, accessibility, distance from the villages, elevation, slope, and positions of transects within land use and landscape and their geographical locations were recorded. Land use regimes were coded SIF (4), CPA (3), TATR (2), and CRUZ (1) based on a perceived disturbance intensity in ascending order. The elevation of each transect within land use was recorded using GPS. A slope/landscape gradient was recorded using Suunto PM-5/369 PC clinometer. These were grouped and coded as follows: 2°-5° degrees (gentle to undulating-coded 1), 5°-6° degrees (moderate-coded 2), 11°-18° degrees (moderately steep-coded 3), and 19°-30° degrees (steepcoded 4) [42]. The position of each transect was recorded as the bottom, foot-slope, mid-slope, shoulder and top as the landform of the terrain was irregular in both forest areas and coded as 1, 2, 3, 4, and 5, respectively [38]. Distance from the village, which may influence the extent of forest and tree species harvesting by the local community [28], was also estimated in km.

Statistical Analysis
The effectiveness of the sampling effort on species richness was evaluated using a species accumulation curve based on Bootstrap estimator in Primer-E [43]. An effective sampling effort captures ≥80% of the estimated species richness [44]. This was then followed by comparing the condition of forest landscape complexity of the two Forest Reserves. First, SIMPER (Similarity Percentage) analysis of land use regimes was done on the original abundance matrix to identify dominant species of each FR [36]. SIMPER also provides an output on the contribution of a species to intra-group (within Forest Reserves) similarity by taking the average contribution of ith species (Av. Sim), of overall pairs of sample plots within a group (j, k), of a species in the Bray-Curtis similarity formula (Equation 1).
where Sjk (i) represents the similarity between the Jth and kth sample, yij represents the entry in the ith row and J column of the abundance data matrix, that is the abundance for the ith species in the jth sample (i = 1,2…p, j = 1,2…n). To describe and compare the spatial patterns of land use disturbance, change drivers for each FRs were subjected to a correlation test using the draftsman plot routine. This is a routine that provides a Spearman correlation coefficient of pairwise variables of all combinations of variables [43].

The Effect of Land Use Regimes on the Difference of Mean β-Diversity in HMFL
To assess the difference in mean β-diversity between land use regimes within each Forest Reserve, first, the information was put into an abundance-based species-sample matrix. The original abundance-based species-sample matrix was transformed into a presence/absence format followed by preparation of the Jaccard similarity coefficient matrix [43]. The Jaccard similarity coefficient calculates the likelihood of a single species being picked at random from two sites without considering the joint absence [36] as follows: where Sjk (i) represents the similarity between the Jth and kth samples and S is the probability (×100). a, b, and c represent the number of species which are present in both samples, the number of species present in sample J but absent from sample in K, and the number of species present in sample K but absent from sample in J, respectively. First, the Jaccard similarity coefficient matrix was then subjected to non-Metric Multidimensional Scaling Ordination (nNMDS) to visually assess the patterns of mean β-diversity of local assemblage among different land use regimes for HMFLs [43,45]. Points in ordination plot with the same color represent plots within the same land use regime, while the different colors represent different land use regimes. The closeness of points to each other in an ordination plot shows that the degree of similarity of sample points in their local assemblages or lower value for the Jaccard dissimilarity coefficient.
The Jaccard similarity coefficient matrix was subjected to Multivariate analysis of variance in PERMANOVA (permutation-based MANOVA) with 999 permutations to test if there is a statistically significant difference (p ≤ 0.05) in the mean β-diversity of local assemblages (group centroids) among the land use regimes. This was followed by a post-hoc pairwise comparison between land use regimes. To measure the effect size of the difference dissimilarity in mean between pairwise land use regimes, distance from pairwise centroids (Av. DJ) was calculated using a Distance among centroids routine. An Av. DJ takes a value between (0, 1), with the ends of the range representing the extreme possibility; (0) represents identical species assemblages between sample points and (1) represents distinct local assemblages between sample points [36]. Changes in mean β-diversity can happen either due to species richness differences (local extinction/immigration), changes in species composition (replacement of species identity), or both between sample points [46].

Species Richness Difference (Δ ) along the Land Use Gradient
To determine if there was a significant difference in species richness between land use regimes, the Margalef index (d) for overall species richness was calculated using the DIVERSE function in Primer-E 7 [36]. The Margalef index (d) is an indicator for species richness/count (S) and takes into consideration the effect of size (N) and the fact that within a larger number of individuals, more species are expected [43].
A resemblances matrix of d-sample plots was then developed using the Euclidean distance. This was analyzed using PERMANOVA (p ≤ 0.05). When PERMANOVA is used to do a univariate ANOVA, the P-values are obtained by permutation, and therefore it avoids the assumption of normality [36]. A pairwise comparison was then used to compare the species richness difference and the conformity to IDH using PERMANOVA. This was done separately for both HMFLs.
Where significant differences were detected, this was then followed by a Hedge (g) metric calculation to detect the effective size of the richness difference between pairs of land use regimes with the same HMFL. The Hedge (g) metric is a weighted average mean standard difference based on a pooled variance measure [45]. It was calculated as follows: where Xa and Xb refer to the mean of paired samples, and SDpooled refers to the pooled standard deviation. SDpooled was calculated as follows: where na and nb refers to the sampling size of the paired samples, while SDa 2 and SDb 2 refer to the square of the standard deviation of the paired samples.
Since Hedge (g) is a biased estimator of population effective size, we used the commonly used J correction factor to calculate the biased corrected Hedges' g value or g * = gJ.

Change in Species Composition (Identity Replacement) along Land Use Gradient
To determine if there were any significant differences of species composition along the land use gradient, PERMDISP (a test of homogeneity of dispersion) procedure was employed on the Jaccard similarity coefficient matrix (p ≤ 0.05). When PERMDISP is used on the Jaccard similarity coefficient, it provides a test of significance between the sampling points on the identity of species they contain. This was then followed by pairwise comparison on species composition between land use regimes. Also, PERMDISP generates a mean square distance of a site to a group centroid (within group dissimilarity) that can directly be interpreted as the percentage of unshared species within a group when the Jaccard dis/similarity matrix is used [36].

A Local Contribution of the Human Modified Forest Landscape to overall β-Diversity of Forest Reserve
To analyze the impact of human modification of forest landscape on overall structure (variability) of local assemblage, a Jaccard coefficient matrix composed of HMFLs and SIF was organized. The three levels, TATR, CRUZ, and CPA were nested under HMFL (Figure 3). After visually inspecting the patterns of overall mean β-diversity of local assemblage among different HMFLs and SIF using Principal Coordinate Principal Coordinated Analysis (PCA), the Jaccard coefficient matrix was subjected to a PERMANOVA test between HMFLs and SIF for each FR. To analyze the impact of human modification of forest landscape on overall mean species richness (đoverall), a species richness (d) resemblance matrix of each FR that included SIF was prepared using the same procedures as above. This was then subjected to a PERMANOVA test between HMFLs and SIF for each FR.
To analyze the contribution to HMFLs and SIF to overall β-diversity of the FR, the Jaccard coefficient matrix that was composed of HMFL and SIF was used. This was subjected to a pairwise PERMDISP test (test of homogeneity of dispersion) between HMFLs and SIF. This was then followed by a PERMDISP test for all land use regimes of the HMFLs to generate the relative importance of each land use regimes contribution to the overall β-diversity of the HMFLs (Mean distance from overall centroid in %). PERMDISP generates the mean square distance of a site (a sample plot) from an overall data centroid (an ideal local assemblage) that can be interpreted as a Local contribuation to β-diversity (LCBD). The distribution of LCBD values between the HMFL and SIF was then analyzed using the Mann-Whitney U test for significant difference as there was a significant difference in dispersion (unequal variance) between sampling units, and further inspected using a Principal Coordinate Analysis (PCO) diagram. In the PCO diagram, the site with high LCBD is that found far from the multivariate centroid of the graphs. If a site has large LCBD, it either indicates a high conservation value due to unique assemblage or conversely may indicate a degraded site with poor species assemblage that may need a restoration action [47].

The Influence of Land Use and Other Environmental Drivers in overall β-Diversity of Human Modified Forest Reserve
Based on the output of the draftsman plot routine of all change drivers, the multicollinearity test of correlation between all pairs of the change drivers was found to be below the acceptable cut-off threshold (R = 0.95). This was followed by Distance-based linear modeling (DISTLM) and distancebased redundancy analysis (dbRDA) to explore the link between β-diversity with change drivers. DISTLM relies on multiple regression models that can accommodate a mixture of categorical and continuous predictors using dbRDA. This is a constrained ordination of sample sites using the same resemblance matrix on the Jaccard similarity coefficient [36].
The relationship between each environmental variable and overall β-diversity of the FRs was initially analyzed separately (excluding other variables) in the marginal test. Variables were then subjected to a forward selection procedure (sequential test, R 2 selection criterion), in which the amount of variability explained by each variable added to the model was conditional of the variables already in the model. P-values for the marginal tests were obtained by using 999 permutations and using the Jaccard similarity coefficient matrix. Distance-Based redundancy analysis was used to visualize the results of the DISTLM [36].

Description of Forest Complexity Condition
The study recorded 2125 number of individual trees in total; out of which 957 and 1168 number of trees were from MFR and TVFR, respectively. The total number of species observed in the whole study area, in MFR and in TVFR were 110, 72, and 88, respectively. The species accumulation curve based on the Bootstrap technique estimated the whole study area, MFR, and TVFR to host about 125, 82, and 99 species, respectively ( Figure 4). Hence, the sampling technique used in this study captured 88.70% of the total species estimated for the whole study, 88.00% for MFR, and 88.18% of TVFR. In MFR, both the observed and estimated species number showed SIF had relatively the highest number of species, followed by CRUZ, CPA, and TATR in descending order (Table 1). In TVFR, both the observed and estimated species number showed that CRUZ had the highest relative number of species followed by TATR, SIF, and CPA in descending order. Similar to the whole study and Forest Reserves, the sampling technique used in this study also captured the majority of the species in all  The SIMPER analysis revealed that the MFR was dominated by seven tree species that contributed about 70% of the total abundance of tree species for the whole landscape, out of which Englerophytum maglismontanum, Bridelia micrantha, and Psidium guajava accounted for about 50% of the total abundance of the trees. In TVFR, 12 tree species dominated and contributed about 70% of the total abundance, out of which Syzygium gerrardii, Xymalos monospora, Englerophytum maglismontanum, Aphloia theiformis, Podocarpus falcatus, and Cassine eucleiformis accounted for 50% of the total abundance of the tree of the landscape. Table 2 shows that there was a negatively strong correlation between land use intensity gradient and distance from the village in both Forest Reserves. The correlation of land use regime with access, the position of the terrain and slope gradients were weak. However, the two Forest Reserves had substantive differences in strength of correlation coefficient with elevation despite having almost a similar mountain range. While land use intensity gradient had a negative and very strong correlation with elevation in MFR, this correlation was weak in TVFR.

The Effect of Land Use Gradient on of Mean β-Diversity of HMFL
The visual inspection of NMDS in MFR ( Figure 5) showed that the local assemblage of all land use regimes in HMFLs was distinct from each other. The PERMANOVA test results in Table 3 show that there was a significant dissimilarity in mean β-diversity among different land use regimes in MFR (F2 = 7.39; P = 0.001). Both pairwise comparison and distance between pairwise centroid (Av. Dj) confirmed that all land use regimes contain highly distinct local assemblage from each other ( Table  4).  Land use gradient (La, fixed factor, three levels) and Transect (Tr, random factor) were nested in La. Degree of Freedom (df), Sum of square (SS), F ratio (Pseudo-F), Permuted probability values (P) are shown. The visual inspection of NMDS in TVFR showed that the local assemblage of all land use regimes is not as distinct as in MFR. Although there was a significant dissimilarity in mean β-diversity among land use regimes in TVFR; (F1 = 3.81; P=0.001), the pairwise comparison in TVFR revealed that there was a significant difference among all pairs of land use regimes. However, the distance between pairwise centroid (Av. Dj) showed that (TATR & CRUZ) was fairly similar. Also, the local assemblages between (TATR & CPA) was at the mid-point of the similarity-dissimilarity continuum (Av. Dj = 0.5).

Change in Species Composition (Identity Replacement) along Land Use Gradients
In MFR, the PERMDISP test showed that there was a significant difference (F2, 36 = 9.11, P = 0.001) in species composition along the land use gradient. The pairwise comparison of PERMDISP result showed that there was a significant difference in the species composition between (TATR, CRUZ) and between (TATR, CPA) with (t = 0.002; p = 0.002) and (t = 2.44; p = 0.003), respectively. The local assemblage in TATR, CRUZ, and CPA had about 37.2%, 49.4%, and 44.9% within group dissimilarity (percentage of unshared species), respectively.

The Local Contribution of the Human Modified Forest Landscape to Overall β-Diversity of Forest Reserve
The PERMANOVA test result showed that there was no significant difference between the mean β-diversity between (HMFLs, SIF) in both Forest Reserves; (F1 = 0.68, P = 0.845) in MFR and (F1 = 0.26, P = 0.817) in TVFR.
The overall mean species richness (đoverall) of the whole MFR was 3.00, and the standard deviation (SD) was 0.46. The mean đ (SD) for HMFLs and SIF were 2.72 (0.64) and 2.82 (0.67) respectively. PERMANOVA test for species richness difference (Δ ̅ ) between HMFL and SIF confirmed that there was no significant difference (F1 = 1.98, P = 2.86). The PERMDISP test revealed that there was a significant difference (F1, 52 = 16.58, P = 0.001) in within group dissimilarity in species composition between HMLF (Mean = 56.6%) and SIF (Mean = 50.5%). Mann-Whitney U test showed the mean rank of LCBD in HMFL (32.67) was significantly higher than the mean rank of LCBD of sample plots in SIF (14.6) (U = 91, P = 0.00). The PCO confirms most of the HMFLs contains many sample plots that were further from the ideal local assemblage of the FR than SIF in FR ( Figure 6). The overall mean species richness (đoverall) of the whole TVFR was 3.75, and the standard deviation (SD) was 0.93. The mean đ (SD) for HMFLs and SIF part were 3.78 (1.02) and 3.66 (0.61), respectively. PERMANOVA test for species richness difference (Δ ̅ ) between HMFL and SIF found that there was no significant difference (F1 = 0.014, P = 0.741). The PERMDISP test revealed that there was no significant difference in within group dissimilarity in species composition between HMLF (Mean = 56.6%) and SIF (Mean = 50.5%) (F1, 57 = 16.58, P=0.331). The Mann-Whitney U test showed the difference in mean rank of LCBD between HMFL (31.16) and SIF (26.6) was not significantly higher than the mean rank of LCBD sample plots in SIF (U = 279, P = 0.37). The PCO shows sample plots in HMFLs and SIF have fairly similar distribution of sample plots from the ideal local assemblage of the FR (Figure 6).

Change Drivers Influencing Overall β-Diversity of Forest Reserves
The marginal test using DISTLM showed that each element of change drivers was found to be statistically significant in explaining the overall β-diversity of both Forest Reserves. The total sum of the individual contribution of each change driver explains about 53.56% and 38.79% of overall βdiversity in MFR and TVFR, respectively. However, the contribution and the significance of those elements in total explain 39.55% in MFR and 28.63% in TVFR of the overall β-diversity when tested with the sequential test of DISTLM (Figure 7). The drop from the total marginal contribution of each change driver in the sequential test may indicate the prevalence of covariance or interaction of land use and other change drivers in shaping the overall β-diversity patterns of a landscape. In MFR, land use regimes explained the highest portion of (13.85%) of the overall β-diversity pattern. This was then followed by elevation (8.42%), position (9.17%), and distance (8.40%). In TVFR, elevation explained the highest variability (9.8%) followed by a position (5.3%) and slope gradient (3.8%).

Discussion
Overall, our sampling technique sufficiently captured the majority of species estimated at all study areas, Forest Reserves, and land use regimes ( Figure 2; Table 1). It is evident that anthropogenic disturbance alters the species diversity through a combination of human land use factors by (i) directly removing preferred tree species for livelihood (e.g., fuelwood and timber); (ii) arresting local successional recovery through the recurrent use, including grazing and herbal medicines; (iii) directly or indirectly changing the local conditions (e.g., soil, moisture, sunlight, competition) [48] and; iv) applying different social norms to govern different parts of a forest landscape [34,35]. Notwithstanding this, the impacts of land use disturbance on biodiversity are neither temporary nor fully avoidable [49]. Furthermore, each landscape, on which different land use intensity is superimposed, is a unique socio-ecological system [26] and differs in their moderating effect on the resilience of local species diversity [30].
The findings from SIMPER and correlational test of land use intensity gradient with other environmental change drivers imply that Mafhela Forest Reserve (MFR) was a relatively simple forest landscape. It had lesser species diversity and dominated by few disturbance tolerant species; on which similar land use regime was spatially clustered across the same elevation range. In contrast, TVFR was complex, species-rich, and dominated by a mix of intermediate and late successional tree species. Unlike MFR, the spatial pattern of land use regimes was more heterogeneous. Hence, a clear understanding of the effect of land use under different conditions of forest landscape complexity in human modified landscapes is crucial for better conservation and management of biodiversity.
Confirming to the hypothesis1, the two HMFLs demonstrated the contrasting effects of land use disturbance to mean β-diversity and conformity to the expectation of IDH. Similar findings of the possibility of divergence and convergence of local assemblages (e.g., reference [21]) and conditional conformity on the response of species richness to intermediate disturbance hypothesis have been reported by many studies (e.g., reference [19]). This can be attributed to the difference in the resilience capacity of the local assemblage of land use regimes due to the variability of landscape moderation to the response of species diversity at different spatial scale [30].
As expected in a simple forest landscape, land use disturbance in MFR enhanced higher dissimilarity in mean β-diversity of local assemblage among all land use regimes. Regardless of the land use intensity gradient, all land use regimes had a distinct local assemblage from each other. This could be because each land use intensity has been creating an environment that suits a set of co-existing species with a similar life-history trait (habitat specialist) [29]. This is consistent with the findings on the lack of conformity of species richness to IDH and the strong evidence of species replacement in HMFLs of MFR. The weak evidence in species richness gradient implies that species replacement played a dominant role in the patterns of β-diversity. Hillebrand et al. [18] underline that when immigration and local extinction (change in species composition) become frequent in a landscape, the species richness can still recover over time to the same level despite the change of species composition.
Furthermore, the higher dissimilarity along a land use gradient in MFR also hints at the limited influence of landscape species pool and the breakdown of ecological connectivity of the landscape (ecological fragmentation) due to the homogenous spatial patterns of land use intensity across the landscape [50]. Ecological fragmentation hinders the recovery of species composition of a vulnerable landscape through natural succession. Conversely, it may enhance further divergence into a set of alternative stable states in the landscape over time [21].
Contrastingly, in a relatively complex forest landscape of TVFR, the local assemblage of land use regimes shared the majority of species among each other. However, the similarity in local assemblage declined along the land use gradient. Such a pattern hints at the higher influence of landscape species pool, ecological connectivity [51], and the positive influence of forest landscape complexity during disturbance-recovery [30]. The conformity of species richness response to IDH in the presence of clear gradient in mean β-diversity hints at the fact that species replacement was practically insignificant to override an orderly local extinction/gain gradient of the TVFR. Hence, the observed β-diversity patterns in HMFLs can be explained by how disturbance affects the mechanisms of species coexistence in a relatively complex landscape.
As expected in the IDH, the intermediately disturbed TATR had a maximum species richness by delaying the competitive exclusion or promotion of co-existence between different life-history traits. The similarity of species shared by CPA and CRUZ with TATR (pairwise DJ), indicates that TATR retained the majority of competitively inferior species of CRUZ as well as the competitively dominant species of CPA that could have become locally extinct under too rare or severe disturbance regimes [7,19]. The severe decline in species richness in local assemblages between (CPA, TATR) may imply that the better traditional protection of CPA resulted in a very severe local extinction of competitively inferior species due to a competitive dominance of a few late successional species.
The findings in TVFR are not surprising considering the recent claim by Munyati and Sinthumule [52] on the decline of deforestation rate and recovery of forest conditions in TVFR vegetation. The conformity of species richness response to IDH may hint that the local assemblage of land use regimes (along a land use intensity gradient) in HMFL in TVFR are more resilient to land use disturbance. However, resilience does not mean the absence of dynamism. Even in the absence of human disturbance, the local neighborhood effect, together with biotic and abiotic elements, may still incur small scale changes in species composition [23].
Overall, the higher local contribution of HMFLs to overall β-diversity in MFR through different land use activities might be contributed substantively to γ-diversity (overall biodiversity) of the landscape. The contribution might also appear as a confirmation for the recent criticism on the insufficiency of some protected areas to cover the scale of compositional dissimilarity (e.g., references [2,13]). However, higher overall β-diversity does not automatically imply that human modification of forest landscape enhances the quality and amount of biodiversity [24]. For instance, the highest contribution of CRUZ in MFR to LCBD implies that the substantial proportion of the HMFL was a degraded ecosystem wherein the substantive parts of the original biota were replaced by competitively inferior and early successional species at the local level.
The fact that land use regimes followed by distance and elevation explain the most substantial proportion of overall β-diversity (Section 3.5) highlights that the current condition of forest complexity at a landscape level is highly simplified by anthropogenic disturbance. As such, the landscape species pool may not prevent the local extinction of old-growth tree species [32] unless it is restored. Instead, in line with the intermediate landscape complexity hypothesis [30], the stateprotected indigenous forests (SIF) appeared to be more effective in safeguarding species-rich hotspots of a vulnerable landscape. However, when considering the findings of Laurance et al. [30] on the impact of environmental deterioration outside of the majority of tropical protected areas for the ecological health of the interior part of protected areas, even the sustainability of the remaining relatively species-rich SIF is uncertain. This is mainly due to the reduction of the minimum dynamic area required for proper ecosystem functioning [52] and a desperate demand for forest products for rural livelihoods. The shift in local species composition may have a detrimental effect on ecosystem provision to the local community.
In a contrasting ecological manner, human modification in TVFR did not adversely impact the overall landscape assemblage and had equal mean local species richness and LCBD. Following the landscape moderated insurance hypothesis [30], the presence of better conditions of forest landscape complexity might have been assisting in a rapid recovery of ecological process under a continuous land use pressure [32]. The fact that natural factors, such as elevation, slope gradient, and position of the terrain are leading drivers to explain the overall β-diversity shows that human modification did not override the natural gradient of γ-diversity. Instead it kept the species diversity in a dynamic equilibrium. Natural factors become ecologically important when the anthropogenic impact disturbance is minor or when disturbances are widespread and fairly uniformly distributed across a landscape [31]. As such, the effectiveness of SIF in TVFR is lower than its counterpart in MFR. The cumulative impact of land use on HMFL did result in comparable average local species richness and a contribution to overall β-diversity with their counterpart SIF.

Conclusions
With the recent prediction of mass extinction of species and the decline of ecosystem services, the debate on whether to maximize on the potential of human modified forest landscapes as an alternative or complementary strategy to protected areas is a non-trivial issue. However, the response for the effectiveness of both conservation alternatives lies in our understanding and response to the question: under what conditions does anthropogenic disturbance enhance, erode, or remain harmless to the pre-existing natural forest conditions in the human biosphere? The contrasting findings between the simple and the complex human modified forest landscape reflect the contrasting insurance value of existing conditions of forest landscape complexity between the HMFLs of the two Forest Reserves. They also reflect that the resilience capacity of local assemblages due to land use pressure can neither be fully explained by one theory nor captured by one species diversity matrix (e.g., Richness or beta-diversity). Using the species richness index alone may obscure the effect of area by averaging the local species richness without discerning the colonization and extinction dynamics. Hence, this situation demands the application of complementary theoretical frameworks and multilevel modeling.
The higher dissimilarity of local assemblages, in the absence of conformity of species richness response to IDH, of a simplified MFR may imply that the replacement of original forest biota in HMFLs is a function of the local extinction of intermediate and late successional species in the whole Forest Reserve. Conversely, a higher species richness of state-protected indigenous forests implies that strictly protected areas can be an effective conservation tool to protect biodiversity hotspots in a simplified forest landscape. They can serve as refugia and source of propagules for the recovery of the local lost species if complemented with restoration efforts for the overall forest landscape complexity. Moreover, reconfiguring the spatial patterns of land use regimes across the HMFLs to increase landscape connectivity may also play a crucial role in restoration. In contrast, maintaining collaborative and holistic landscape management in TVFR, using IDH as a guiding tool, may ensure the sustainability of the current forest landscape complexity and the retention of the rich species diversity. However, conclusive remarks cannot be made based on beta-diversity, since it is only exclusively dependent on presence/absence data. Further research is needed using abundance based data to explictly expose which species traits and their relative abundance would be most affected by any disturbance. This will help to efficiently allocate increasingly limited conservation resources for the conservation of priority species and habitats. Author Contributions: M.G.A. designed the research, collected data and analyzed under the suprvision of P.W.C. and E.S.P.A. The first draft manscript was written by M.G.A. and P.W.C. and E.S.P.A made a substantive contribuation in the revison. All authors have read and agreed to the published version of the manuscript.

Funding:
The research was funded by the Forest Programme at the University of Pretoria in South Africa through the South African Forest Company Limited (SAFCOL) Forest Chair.