The CBNP is situated on Cat Ba Island, in the Tonkin Gulf, northern Vietnam. It lies ca. 45 km east of Haiphong City and 150 km south-east of Hanoi. Geographical coordinates are 20°44'N to 20°55'N (DM), 106°54'E to 107°10'E (DM) (Fig. 1) (Pham et al. 2020).

Figure 1.

Location of Cat Ba National Park (CBNP), Haiphong, Vietnam. We collected data at three different study sites: MSA, a mid-slope area; LLA, a lowland area; ISA, an isolated area; 4, 5 and 6 are other restricted areas in the CBNP that were not considered for this study (Pham et al. 2020). Map data copyrighted by OpenStreetMap contributors and available from https://www.openstreetmap.org (CC BY-SA 2.0).

The CBNP is located on a limestone archipelago that consists of 366 islands and islets (CBNP 2005). The total CBNP area is 16,197 ha with ca. 10,932 ha of terrestrial island areas and 5,265 ha of marine areas (CBNP 2005, 2007). The average elevation is between 100 and 150 m a.s.l. with the highest peak of 331 m a.s.l. and average slopes ranging from 15° to 35° (CBNP 2007). The CBNP is characterised by two distinct seasons: the rainy season (from May to October) and the dry season (from November to April) with yearly rainfall from 1,500 to 2,000 mm. Average annual temperature is 23 °C and the average air humidity is approximately 86% (CBNP 2005; Le 2006).

Various ecosystems and forest types can be found in the CBNP and these include evergreen forests on limestone, wetland forest ecosystems in steep mountain valleys, mangrove forests, coral reefs and cave systems (CBNP 2005; Le 2006). The forest ecosystems include primary (undisturbed by direct human activity) evergreen broadleaf tropical rain forests, secondary evergreen broadleaf tropical forests in the lowlands (previously disturbed by human activity) and on limestone mountains, secondary moist evergreen restoration forests on limestone mountains and valleys, restored bamboo forests, wetland forests in the limestone valleys, mangrove forests and plantation forests (Carle and Holmgren 2003; CBNP 2007).

In 2004, Cat Ba Island was designated a UNESCO biosphere reserve because of its diverse flora and fauna (CBNP 2005). In total, 1,561 vascular plant species belonging to 842 genera and 186 families have been recorded across the different terrestrial ecosystems. Amongst them, 408 tree species have been observed in the CBNP in the past years by local authorities (CBNP 2007). Twenty-nine tree species listed in the IUCN Red List were recorded in the CBNP (Le 2006; CBNP 2007). Forty-three additional tree species are listed on the Vietnam Red List in terms of conservation (Le 2006; CBNP 2007). To allow the natural development of ecosystems and to protect precious plants and animals in the CBNP, six strictly protected zones were established, totalling 4,915 ha (Fig. 1). In these protected areas, management activities, such as timber logging, hunting or slash and burn, are prohibited. We chose three of the six restricted areas for our study that represent the three main conservation areas of the National Park. These three areas are characterised by primary and secondary evergreen broadleaf tropical forest ecosystems (Fig. 1). These areas also represent an accessibility gradient. The lowland area (LLA) (1,916.4 ha) is in the centre of the Park and is the most easily accessible of the three. The mid-slope area (MSA) (600 ha) lies in the north-western part of the Park and is characterised by steeper slopes compared to the LLA. The third area is an isolated area (ISA) (1,557.8 ha) and is located on a separate island in the eastern part of the Park. Its island situation isolates ISA from the two other study areas (Fig. 1). Note that in Fig. 1, some local villages were included into the boundaries of MSA, making the study area larger than the actual forested area.

We used the simple random sampling technique (Kleinn et al. 2009) to set up sample plots (Fig. 2). Thirty strips were created in each research area. In each strip, random sample plots were generated using random numbers to determine their coordinates. Two uniform random numbers U_1i, U_2i (the U interval from 0 to 1) were used each time to calculate X_i = U_1i * X_max, with Y_i = U_2i * Y_max as coordinates for each random sample plot, where X_max and Y_max were the largest coordinates of the area map (Fig. 2). If the coordinate (X_i, Y_i) appeared within the defined strip, this point was accepted as a sample plot point. Otherwise, the point was rejected and the procedure was repeated with two new U(s) random values (Fig. 2).

Figure 2.

Sampling design scheme.

Based on this technique, a total of 90 random sample plots were created (30 plots in each of the three protected areas (LLA, MSA and ISA)). Each sample plot area was 500 m² in size (20 m × 25 m).

To analyse plant species composition in the CBNP and to identify different forest communities, we used hierarchical cluster analysis with Ward’s method using the function ‘vegdist()’ to create a Bray Curtis distance matrix, the function ‘hclust()’ to conduct a cluster analysis and the function ‘cutree()’ to draw a community dendrogram. These functions are implemented in the ‘vegan’ package in R (Oksanen et al. 2019).

To display the spatial distribution of communities in multidimensional space, we used the non-metric multidimensional scaling (NMDS) ordination method, based on abundance data. In the first step, we used the ‘metaMDS()’ function to run an NMDS and then used the ‘envfit()’ function to add environmental data into the NMDS graph using ‘vegan’ in R . We included all vegetation layers into the NMDS. Tree species of the overstorey and regeneration layers were combined to avoid duplication of species names.

To determine indicator species for identified forest communities, we used the function ‘multipatt()’ in the package ‘indicspecies’ (De Cáceres and Legendre 2009). The indicator species analysis calculates the specificity of a species as the number of occurrences of a species within a community relative to the number of occurrences across communities and the frequency of the species within a community as the relative number of occurrences per community. Multiplying specificity and frequency results in an indicator value between 0 (species not occurring in a specific forest community) and 1 (species occurring only and always in a specific community). The most significant indicator species identified amongst the tree and non-tree species were used to name the forest communities that were determined by cluster analysis.

To determine factors affecting community distribution in the CBNP, we correlated the NDMS axes’ values with the environmental factors and human impact indicators. The function ‘anova()’ following a post-hoc Tukey-test with function ‘glht()’ in ‘multcomp’ package was used to compare these variables amongst communities.

We additionally applied different selection operators to build decision trees to weigh the predictors that characterised plant species communities (= the response). Forward selection of the predictors performed best. This method starts with an empty model and adds predictors to explain the response. Performance is estimated in each round using cross-validation. Only predictors with the highest performance increase are kept, then a new round is started. The maximum number of attributes to add was limited to seven to avoid overfitting; the iteration was aborted when performance no longer increased. This yielded a decision tree classifying the plant species communities. Classification trees represent a robust, non-parametric, binary procedure that partitions variance in the variation of communities through a series of splits in more homogeneous groups based on environmental factors or human impact variables (De’ath and Fabricius 2000; Cutler et al. 2007).

We contrasted plant species richness amongst the detected forest communities, species groups (tree and non-tree-species) and vegetation layers (overstorey tree, regeneration, non-tree layer) at different spatial scales. For alpha diversity, we compared plot-based species richness amongst communities using the ‘anova()’ function, followed by a post-hoc Tukey-test with function ‘glht()’ in the ‘multcomp’ package. To compare total species diversity in the identified communities, we used the ‘iNEXT’ package in R (Hsieh et al. 2016) to estimate the respective gamma diversity. The ‘iNEXT’ package applies different orders of Hill numbers (q) and quantifies diversity by applying rarefaction and extrapolation methods (Chao et al. 2014). We considered the first three Hill numbers referring to species richness (q = 0), the true diversity of the Shannon Index, which is the exponential of the Shannon Index (q = 1) and Simpson diversity (q = 2) (Chao and Jost 2012; Hsieh et al. 2016). We estimated diversity for each defined forest community separately and across communities as an estimate for the gamma diversity of the CBNP.

To investigate how much each detected forest community contributed to the gamma diversity of the CBNP, we contrasted the estimated species pool of different combinations of communities to the species pool estimated for a combination of all detected communities. We also estimated the gamma diversity of tree and non-tree species of the CBNP and of different forest communities.

The relationship between non-tree species and the overstorey tree layer and between non-tree species and the regeneration layer was checked using linear regression models. We investigated the relationships across communities and for each community separately.

All statistical analyses were conducted using the statistical software R version 3.4.2 (R Core Team 2017). The level of significance for all statistical inferences was defined as p-value < 0.05.

We identified four main communities that differ distinctly in species composition in Cat Ba National Park, (Fig. 3, Appendix 1). Each community has characteristic indicator species (Table 1). The most reliable indicators (highest indicator index value) amongst the tree and non-tree species provided the community names (Table 1, species characteristics are provided in Appendix 2). We defined the Saraca dives + Calamus tetradactylus community (SCt), the Sterculia lanceolata + Chloris barbata community (SCb), the Ficus superba + Acanthus ebracteatus community (FAe) and the Clausena excavata + Desmos cochinchinensis community (CDc) (Appendix 2). Most indicator species (14 species) were associated with the SCt community, six species were associated with the SCb community and only four indicator species with the FAe community. For the CDc community, ten indicator species were identified (Table 1).

Figure 3.

Spatial distribution of communities in NMDS analysis. The ellipses show the confidence levels at 95% of covariance of each community.

Ordination revealed that the SCt community was found mainly in the lowland area (LLA). The SCb community characterised a transition from the lowland to mid-slope area. The FAe community was associated with the mid-slope area and the CDc community with the isolated area (Fig. 3).

Forest structure differed amongst the four communities. On average, we found the lowest diameter at breast height (DBH), height (Ht), basal area (BA) and volume (Vol) in the FAe community, while the SCt community was characterised by the highest mean volume, DBH, basal area and height. The SCb community had the highest tree species richness, whereas the CDc community was the most species-poor in the overstorey layer (Table 2).

Community comparisons, as well as the decision tree analysis, showed that the communities FAe and CDc were associated with steeper slopes and a high percentage of rock surface area (Fig. 4, Table 3, Appendix 3). In contrast, the SCb and SCt communities were found in areas with deeper soils (Fig. 4, Table 3, Appendix 3). Although soils are shallower in the habitats of the FAe and CDc communities, their soil humus content (S_SH), soil moisture (S_SM) and cation exchange capacity (S_CEC) were high (Fig. 4, Table 3, Appendix 3) in the upper soil layers. Light conditions differed little between the four communities. Overall, SCb received the most light (Table 3, Appendix 3). Indicators of human impact were mainly associated with the communities SCt and SCb (Fig. 4).

Figure 4.

Correlation of environmental and human impact factors if significant (T = topographic variables, S = soil variables, H = human impact) with the NMDS axes. Length of the arrow indicates the strength of the correlation. Arrows in the same direction indicate a positive correlation of the variables with each other and a unidirectional response to NMDS axes. All variables that significantly correlated with the axes (p-value < 0.05) are shown. ISA = isolated area; LLA = lowland area; MSA = mid-slope area. The abbreviations for environmental and human disturbance factors are defined in Table 3.

Differences in environmental factors among communities were also reflected in the differences in environmental factors among the study sites in which the communities were mainly located (see Figs 3, 4, Appendices 1, 4). For example, the average slopes in ISA and MSA were steeper than in LLA. As the FAe and CDc communities were associated mainly with the MSA and ISA sites, they were also characterised by higher slope values compared to the other two communities (Appendix 1). In addition, the ISA site had the highest percentage of rock surface as did the CDc community (Appendix 3).

Across communities, the non-tree species layer was, on average, more species-rich than the overstorey and regeneration layers, with the highest non-tree species richness found in the SCb community (Fig. 5). The overstorey tree layer was significantly more diverse in the SCt and SCb communities compared to the FAe and CDc communities (Fig. 5). A similar pattern was found for the regeneration layer (Fig. 5). Thus, the communities with the lowest plot-based species richness of the regeneration layer also had the lowest plot-based species richness in the overstorey tree layer.

Figure 5.

Boxplots of species richness of the three investigated layers in each forest community. Letters indicate statistically significant differences amongst communities. SCt; Saraca dives + Calamus tetradactylus community. SCb; Sterculia lanceolata + Chloris barbata community. FAe; Ficus superba + Acanthus ebracteatus community. CDc; Clausena excavata + Desmos cochinchinensis community.

Gamma diversity estimations identified significantly higher values for the communities SCt and SCb compared to the communities FAe and CDc. For the Hill number q = 2, the SCb community, which characterised the transition from the lowland to the mid-slope area, was significantly more diverse than the other three communities; confidence intervals do not overlap. For all Hill numbers, the forest community associated mainly with the isolated area (CDc) had the lowest gamma diversity (Fig. 6). When plots of all four communities were combined, gamma diversity did not exceed the most diverse community, indicating that community assemblages are not complementary, but exhibit nested diversity.

Figure 6.

Estimated gamma diversity for different Hill numbers and forest communities. q = 0: species richness; q = 1: Shannon diversity; q = 2: Simpson diversity. The red colour represents the estimated gamma diversity when all communities are considered for estimation. Yellow-green: Saraca dives + Calamus tetradactylus community (SCt). Green: Sterculia lanceolata + Chloris barbata community (SCb). Light-blue: Ficus superba + Acanthus ebracteatus community (FAe). Pink: Clausena excavata + Desmos cochinchinensis community (CDc). Graphs were extrapolated to a sample size of 90 plots.

When investigating the importance of each community to CBNP gamma diversity, we first estimated the species pools of different combinations of communities. We found that the SCb and SCt communities were most important for regional gamma diversity. When these communities were eliminated from community combinations, gamma diversity was greatly reduced as compared to the total species pool and to other community combinations (Appendix 5). If, however, the CDc and FAe communities were not included in the estimation procedure, the estimated species pool was not significantly different from the total species pool, indicating that these communities are mainly subsets of the species-richer communities (Appendix 5).

Gamma diversity patterns of the different communities were driven in large part by tree species (Fig. 7). Tree species were significantly less diverse in the communities FAe and CDc compared to the other two communities (Fig. 7, Appendix 6). For the non-tree species, the difference between the two community groups (SCb, SCt vs. FAe, CDc) was not as pronounced for q = 0 (Fig. 7). For q = 1, the SCb community was significantly more diverse than the communities FAe and CDc. For q = 2, SCb was more diverse than all other identified communities (Fig. 7, Appendix 6).

Figure 7.

Gamma diversity patterns of tree species and non-tree species of the different communities. The graphs show species richness (Hill number q = 0, including 95% confidence interval). Red: all four communities, Yellow-green: Saraca dives + Calamus tetradactylus community (SCt), Green: Sterculia lanceolata + Chloris barbata community (SCb), Light-blue: Ficus superba + Acanthus ebracteatus community (FAe), Pink: Clausena excavata + Desmos cochinchinensis community (CDc). For Hill numbers q = 1 and q = 2, see Appendix 6. Graphs were extrapolated to a sample size of 90 plots.

Tree and non-tree species contributed nearly equally to the diversity of the CBNP (Table 4). However, in the SCt and SCb communities, tree species diversity was higher than non-tree species diversity, while in the other two communities, the converse was true (Table 4).

We did not find a significant relationship between species richness of the non-tree layer and that of the tree regeneration layer. We also found no effect of species richness of the overstorey tree layer on species richness of the non-tree species layer or the regeneration layer (Fig. 8).

Figure 8.

Linear regression models of plot-based species richness relating: a species richness of the non-tree species layer and species richness of the regeneration layer b species richness of the overstorey tree layer and species richness of the non-tree species layer and c species richness of the overstorey tree layer and species richness of the regeneration layer. The shaded area indicates the 95% confidence interval.

Across all communities, there was also no significant effect of the coverage of the non-tree species layer on the species richness or abundance of the regeneration layer (Fig. 9a, c). However, at the community level, we found that coverage of the non-tree species layer negatively affected both abundance and richness of the regeneration layer in the SCb community and of richness in the SCt community. In contrast, coverage by the non-tree species layer positively influenced the richness of the regeneration layer in the FAe community (Fig. 9b, d).

Figure 9.

Linear regression models of the coverage of non-tree species and species richness and abundance of the tree regeneration layer. Graphs (a) and (c) show the relationship between the coverage of non-tree species with the abundance and species richness of tree regeneration across all communities, graphs (b) and (d) for the four communities. Dashed lines show statistically non-significant, solid lines show significant correlations. To improve visibility, we did not include the confidence intervals for (b) and (d).

We did not observe any significant effect of species richness or coverage of the non-tree species layer on tree abundance or richness of threatened tree species in the regeneration layer (Appendix 7).

We distinguished four forest communities in the CBNP that differ in forest structure and are characterised by different environmental conditions (Fig. 3, Tables 1, 2, Appendix 1). Topographic variables, in combination with soil conditions, turned out to be important factors influencing the differences in plant species composition across the CBNP. Topographic heterogeneity can result in different microclimatic conditions, affecting temperature and soil moisture at small scales, thus influencing the abundance and distribution of plant species (Vivian-Smith 1997; Seta et al. 2018). In addition, topographic heterogeneity indirectly acts on plant species composition and diversity through drainage and by differences in the hydraulic regime, which results in different growing conditions for plants on ridges, slopes and in valleys (Trichon 1997). These characteristics may also affect seed erosion and seed accumulation, leading to differences in plant species distributions (Vivian-Smith 1997; Holl et al. 2013).

We found that the SCt and SCb communities occurred mainly on sites characterised by deep soils, whereas the FAe and CDc communities were found on sites with rough terrain (steep slopes and high rock surface), but with relatively high soil nutrient content (high soil humus content and CEC) (Fig. 4, Table 3, Appendices 1, 4). Our environmental characterisation of the communities corresponded well to the investigated study sites. The SCt and SCb communities characterised the lowland area in transition to the mid-slope study area, while the CDc community was associated mainly with the isolated area, characterised by steep slopes. Our results are in line with Zhang et al. (2021), who identified topography as the main factor determining niche differentiation of tree species in heterogeneous tropical limestone forests of China. Topography influenced soil depth and, thereby, most likely determined water availability and drought stress, factors that may have shaped species distribution.

Topographic heterogeneity can also influence light availability, an important determinant of plant species distribution (Zhu et al. 2016). It also affects forest composition by influencing regeneration patterns (Tateno and Takeda 2003). However, in our case, the indirect site factor (L_ISF), which we used as a proxy for light availability, did not affect plant species distributions (Fig. 4, Table 3, Appendix 4). It is important to mention that the gradient of light availability recorded across our study sites was quite low, indicating a relatively homogeneous canopy cover across the CBNP.

Although indicators of human impact could not be clearly connected to plant species distribution, observed differences in the forest communities of the different study areas could have reflected human influence, which was strongest in the lowland area. In the lowland area, with the communities SCt and SCb, species richness was highest (Figs 4, 5, Table 3). In contrast, other studies found negative effects of human disturbance on species diversity; humans extract certain tree species with subsequent effects on species composition (Blanc et al. 2000; Do et al. 2017). Chazdon (2003), however, concluded that forest recovery could be rapid even after large-scale human disturbances when soil and aboveground vegetation were not heavily impacted. In addition, Dao and Hölscher (2018) found a positive correlation between the occurrence of footpaths and the density of tree species used for non-timber forest products. This suggests a positive effect of human activity on some early-successional species. Other studies have confirmed the strong potential of tropical forests to recover after historical disturbances (Lusk and Smith 1998; Bayliss-Smith et al. 2003), but also underscore the role of historical disturbances, whether anthropogenic or natural, in shaping species composition.

We identified several indicator species for the four forest communities that provided valuable information on species-environment relationships. While species of the communities SCt and SCb seemed to be restricted to deep soils, indicator species in the FAe and CDc were able to tolerate harsh soil conditions, such as a high percentage of rock surface and shallow mineral soils (Guo et al. 2017). They were also able to overcome potential dispersal filters due to the island isolation of the ISA. The bird-dispersed Clausena excavata, for example, is characterised by rapid germination, thereby avoiding the risk of desiccation in shallow soils. High seedling survival under various environmental conditions makes this species both a successful invader outside its native range and a successful coloniser (Vieira et al. 2010). Supporting results by Santo‐Silva et al. (2021) demonstrated that isolation favours the abundance of disturbance-adapted pioneer tree species. The annual herb Blumea lacineata, another indicator in the CDc community, is also characterised by an effective dispersal strategy and is classified as a weed species (Wester 1992). On the other hand, Sterculia lanceolata, an indicator species of SCb, can be characterised as a mid- to late-successional species on deeper soils (Zhang et al. 2013). Its co-existence with the tropical weed species Chloris barbata (Holm et al. 1979) in the CBNP suggests that this community and its assemblage of indicator species characterise a successional stage after disturbance in lowland and mid-slope areas. Lower DBH and tree heights in this community, as compared to the SCt community of the lowlands, support this assumption. The species assemblage of the Saraca dives + Calamus tetradactylus reliably characterises lowland sites, as Saraca dives is a dominant species in foothills and in valleys with high and constant water availability in limestone forests (Zhang et al. 2021).

Thus, our study underscores the value of using identified indicator species or groups of species for an overall assessment of the environmental conditions in tropical limestone forests of Southeast Asia. By monitoring certain indicator species, shifts in environmental conditions can be reliably detected.

We identified in total 302 species belonging to 112 families in the CBNP (Table 4, Appendices 5, 8) within the 90 recorded sample plots. The total species pool was estimated at 368 species. We found that, as compared to other studies, plant species diversity in the CBNP was quite high. For example, diversity was higher than that of transitional rainforest vegetation in southwestern Ethiopia (from 130 to 139 species) (Assefa et al. 2013), neotropical primary and secondary forests (100 species found on two 1 ha plots (van Andel 2001), or subtropical forests in China (240 species in 27 plots of up to 900 m² (Both et al. 2011). It was less diverse, however, than rainforests in Columbia (442 vascular plant species in 0.9 ha; (Galeano et al. 1998) or seasonal dry tropical forests in Mexico (917 species in 20 representative sites of 0.1 ha; (Trejo and Dirzo 2002). Numbers are comparable with forests of Dinagat Island in the Philippines where 432 native plant species have been recorded (Lillo et al. 2019). In Malaysian tropical rain forests, 825 species in 50 ha forest inventory plots have been detected (Fangliang et al. 1997). Although differences in species numbers among the aforementioned studies may also be due to different experimental designs and scales of inventory (Ferraz et al. 2004; Cicuzza et al. 2013; Júnior et al. 2014), our results indicate that plant species diversity in the CBNP is comparable to other tropical and subtropical forests worldwide. Interestingly, non-tree species contributed roughly 50% to total species diversity (Table 4) emphasizing both the importance of different growth forms for biodiversity and confirming other studies (Gentry and Dodson 1987; Linares-Palomino et al. 2009). The forest communities FAe and CDc, characterizing rough and rocky terrain, had a higher share of non-tree species than those of the lowland and mid-slope areas, indicating that this species group is shaped by stochastic assembly processes, while the distribution of tree species is mainly driven by environmental differences among habitats (Both et al. 2011). Similarly, isolation seems to have acted as a stronger barrier for tree than for non-tree species (Hill and Curran 2003; Martin-Queller et al. 2017).

In contrast to our expectations, the heterogeneity we found in forest community composition did not appear to drive overall plant species diversity of the CBNP. The linkage between isolation and species diversity has been investigated in previous studies (Quinn and Harrison 1988; van Andel 2001; Slik et al. 2003; Scheffer et al. 2006; Chytrý et al. 2010). Isolated conditions can promote species diversity and the distinctiveness of communities (Ferraz et al. 2004), and can also drive speciation, leading to species endemism (EI-Bana 2009; Chytrý et al. 2010). Our findings indicate that the isolated area has the lowest species diversity in our study (Figs 6, 7, Appendix 5). This supports the conclusion by van Andel (2001) that relative isolation can also have a negative effect on species diversity. Slik et al. (2003) confirmed that diversity was negatively correlated with isolation and habitat size in North-eastern Borneo forests. Our results, therefore, suggest that rough terrain conditions in combination with the isolated location function as a strong dispersal and establishment filter that only a subset of the plant species can overcome. Such species are characterized by good dispersal ability (see above). The missing complementarity in species diversity among communities, however, indicates that the species assemblages colonizing the isolated area are instead a subset of richer species communities. Species colonizing the isolated sites can therefore be characterized as generalists growing under various site conditions (Vieira et al. 2010; Santo‐Silva et al. 2021).

Our results also confirm the theory that larger areas boost high species diversity. The two species-rich SCt and SCb communities (Figs 3, 6, Table 4, Appendix 5) were located mostly on the main island of the CBNP; this presumably provided more habitats and better habitat connectivity with positive effects on the dispersal ability of species (MacArthur and Wilson 1967; Hill and Curran 2003; Martin-Queller et al. 2017). This agrees with Scheffer et al. (2006), who stated that local diversity is reduced in response to isolation through dispersal limitation. They also implied that larger local habitat patches and greater connectivity facilitate regional biodiversity. Turner et al. (1996) even found that a long period of isolation in lowland tropical forests in Singapore was a major contributor to species loss.

The species-rich SCb and SCt communities were associated with greater soil depth, while the percentage of rock surface and slope were correlated with species-poorer communities (FAe and CDc communities) (Fig. 4, Table 3). Slope and rock outcrops have a huge effect on soil cover in our research areas (Chytrý et al. 2010; Seta et al. 2018), influencing soil nutrients (Chytrý et al. 2010) and water availability (Zhang et al. 2021). This seems not only to result in differences in community composition as discussed before, but also in different numbers of tree and non-tree species. Tree species diversity was more strongly influenced than non-tree species by differences in topography, with more tree species colonizing the lowlands. The Simpson diversity of non-tree species was particularly high in the SCb community. Here the presence of indicator Chloris barbata suggested an impact of former disturbance. We do not have, however, enough information on former anthropogenic or natural disturbance events to make a definitive connection.

Humans can have a huge influence on floristic diversity; there are numerous examples from around the world. For example, Júnior et al. (2014) found that human activities changed tree species diversity and forest structure in semi-deciduous forests in Brazil. Assefa et al. (2013) found that anthropogenic interferences reduced species diversity in southwest Ethiopian forests. In our case, indicators of human impact and species-rich forest communities were positively associated (Fig. 5, Table 2). These results may indicate that the natural recovery potential of a species after historical disturbances is high (van Andel 2001; EI-Bana 2009). In any case, reconstruction of former disturbances would be useful to evaluate the potential for tree species recovery as well as to understand the effects of former disturbances on current species composition. However, as differences in human impact and environmental conditions among communities and study sites interact, we cannot verify a potential impact of former human disturbance on plant species diversity and composition in the CBNP.

Across communities, we found no interactions amongst forest layers (Both et al. 2011), but found contrasting patterns for the different communities. In the two species-rich communities (SCb and SCt) found on deep soils, we detected negative effects of non-tree species coverage on tree species richness and abundance in the regeneration layer (Figs 8, 9b). It is well known that tree regeneration competes with herbaceous species on sites with sufficient water and nutrient availability (Both et al. 2011). Interestingly, under harsh environmental conditions in the CDc communities and also in the FAe (Fig. 9d), coverage by non-tree species appeared to have neutral (CDc) to positive (FAe) effects on the tree regeneration species richness; this supports the hypothesis that changes in water availability, such as those induced by more shallow soils, can shift plant interactions from competition to facilitation (Bertness and Callaway 1994; Holmgren et al. 1997). Along the steep slopes with high rock coverage, non-tree species can improve conditions for tree regeneration, for example, by trapping litter and humus. This, in turn, can improve water and nutrient status, supporting seedling germination and survival (Yirdaw et al. 2015; Yirdaw et al. 2019). This may also explain higher soil moisture values found in the topsoils of the FAe and CDc communities as compared to SCt and SCb. In addition, litter and humus accumulation may generally improve seed storage, preventing rapid erosion.

Competition between tree regeneration and non-tree species in the SCt and SCb communities may also have resulted from former canopy disturbances. Increased light levels induce understorey growth with potential negative effects on tree regeneration. The light-demanding weed species Chloris barbata, identified as an indicator species in the SCb community, may be an indication of former disturbances as are the higher values of human impact indicators in SCt and SCb. However, the light availability measured for this study was around 9.17% ± 6.4% (mean ± sd) and homogeneous across communities.

Although we have considered many environmental factors and indicators of human impact to explain plant species composition and diversity in the CBNP, we acknowledge some shortcomings. We did not explicitly consider climatic factors (Slik et al. 2003; Cicuzza et al. 2013), such as rainfall and temperature. However, since sites were close to one another (compare Fig. 1), we did not expect great differences in climatic conditions. We also do not have detailed information on the disturbance history of the CBNP, even though evolutionary or historical disturbance processes may control the distribution of certain species in certain sites (Trejo and Dirzo 2002; Tuomisto et al. 2003; Chytrý et al. 2010). Finally, factors such as geographical distance or tectonics may also have played a role in the observed patterns of plant species composition and diversity; both of these factors might affect environmental conditions, habitat formation and floral assemblages (van Andel 2001; Slik et al. 2003; Júnior et al. 2014). In future research, it may also be worth focusing on plant functional groups, which would require more research into the autecology of the different tree and non-tree species.