Three parameters, namely the uniform angle index, dominance, and mingling (Li et al., 2012) (Table 1), were selected to analyze the stand spatial structure of the northern tropical karst seasonal rainforest from 2011 to 2021. Spatial structure was characterized using zero, univariate, bivariate, and trivariate distribution of these parameters. The zero distribution represents the scenario in which each of the three spatial structural parameters take its mean value. The univariate distribution describes the frequency of five possible value levels for any single structural parameter. The bivariate distribution was derived by systematically pairing structural parameters across their ordinal levels, generating 25 distinct relative frequency distributions. The trivariate distribution involves the cross-classification of five value levels for two parameters (X and Y) with the five value levels of a third parameter (Z), yielding 5×5×5=125 spatial structural combinations and their corresponding relative frequency distributions.

Relative frequency was calculated by dividing the number of individuals in each value or value combination by the total number of individuals. The uniform angle index, dominance index, and mingling index were each classified into five distinct value levels: 0, 0.25, 0.5, 0.75, and 1. For the uniform angle index, these values correspond to very uniform, uniform, random, irregular, and very irregular distributions, respectively. For the dominance index, the five values represent predominant, subdominant, intermediate, disadvantaged, and absolutely disadvantaged. For the mingling index, the five values correspond to zero mingling, weak mingling, medium mingling, high mingling, and complete mingling. To minimize edge effects, a 5 m buffer zone was applied to the plot.

To quantify individual dynamics, live trees, recruited trees, and dead trees were identified between 2011 to 2021. Sankey diagrams were then employed to visualize transitions among these categories, and the proportions of change in the three processes were calculated to reveal dynamic shifts in stand spatial structure at the individual tree level. Tree recruitment was defined as individuals that either germinated after 2011 or had DBH < 5 cm in 2011 but grew to DBH ≥ 5 cm by 2021. Tree survival was defined as individuals that were alive in both 2011 and 2021. Tree mortality referred to individuals with DBH ≥ 5 cm in 2011 that had died by 2021.

The values of surviving trees exhibited three possible trends: “ascending”, “stable” and “descending” (Zhou et al., 2022). For example, an increase in value from 0 to 0.25 was defined as “ascending one-step (Flow+1)”, whereas a decrease from 1 to 0 was defined as “four-step descending (Flow-4)”. Cases in which the value remained unchanged were defined as “stable”. In total, nine possible value-change scenarios were recognized for surviving trees. By comparison, five value-change scenarios were applied to tree recruitment and mortality processes.

We used the forestSAS (https://CRAN.R-project.org/package=forestSAS) package in R 4.3.1 (https://www.r-project.org) to analyze the stand spatial structure, and the ggalluvial package to draw Sankey diagrams.

Spatial structural parameters and environmental factors were calculated for each 20 m×20 m sample plot. The environmental factors primarily include topographic factors such as elevation, aspect, slope, convexity (Harms et al., 2001; Valencia et al., 2004), topographic wetness, elevation above channel, and rock-bareness rate (RBR). Calculation methods for these indices followed Punchi-Manage et al. (2013) and Kanagaraj et al. (2011). Aspect was transformed using the following formula:

where TRASP denotes the transformed aspect index, and aspect is the slope direction angle (°). This transformation normalizes TRASP values between 0 and 1. Higher TRASP value indicate drier and hotter habitats: a TRASP value of 0 corresponds to the north-northeast direction, whereas a TRASP value of 1 corresponds to the south-southwest direction.

The factors included maximum diameter at breast height, the sum of the area at breast height, the number of individual plants (Ind), and mean diameter at breast height. When calculating structural parameters for each 20 m×20 m plot, a 2 m buffer zone was applied.

To determine whether the stand spatial structure exhibited spatial autocorrelation, we performed a global spatial autocorrelation analysis using the following calculation formula:

where S₀ is the sum of all elements in the spatial weights matrix; n is the number of spatial units, x_i and x_j are the attribute values of unit i and unit j respectively, x̄ is the mean of all attribute values, and w_ij is the spatial weight between units i and unit j (de Campos et al., 2024).

Based on the spatial dependence of stand spatial structure data and the objectives of this study, both the spatial lag model (SLM) and spatial error model (SEM) were constructed to systematically analyze spatial autocorrelation and identify the key driving factors of stand spatial structure (Wang et al., 2019). The SLM is expressed as:

The SEM is expressed as:

where Y is the dependent variable matrix, X is the explanatory variable matrix, W is the spatial weight matrix, β is the vector of regression coefficients, and ϵ is the model residual. The terms ρ, λ, and γ represent spatial autoregressive coefficients, and μ denotes the spatially correlated residual.

In the global spatial autocorrelation analysis, the mingling index and uniform angle index exhibited significant positive spatial correlations, whereas the dominance index showed no significant spatial autocorrelation. Furthermore, results from both the SLM and SEM yielded non-significant, indicating that neither the spatial structure of the stand nor its influencing factors showed significant spatial autocorrelation (Table 2).

Then the impact of environmental factors on the spatial structure of the stand was analyzed using generalized additive models (GAM). The GAMs are well suited for fitting non-linear relationship between the response variable and multiple explanatory variables. In this study, the smoothing function was specified as the default thin-plate spline with degrees of freedom set to 8. Prior to model fitting, pairwise correlation coefficients were calculated among the explanatory variables, and a threshold of |r|>0.7 was used to identify potential multicollinearity. Variance inflation factors (VIF) were then calculated, with values < 10 indicating no substantial collinearity. Model reliability was evaluated using the deviance explained, while the significance of explanatory variables was assessed using F-tests. To further evaluate the relative importance of explanatory variables, the fitted GAM was decomposed using the gam.hp package to obtain contribution rates (Lai et al., 2024). The expression for the GAM is:

where E(Y) is the expected value of the response variable, G[] is the link function determined by the distribution of the response, β₀ is the intercept, and f₁, …, f_n are smooth functions of the n explanatory variables, usually estimated using smoothing splines (Wood, 2006).

In this study, four diversity indices were used to characterize sapling diversity: the Shannon-Wiener index (H), the Simpson index (D) (Simpson, 1949), the Pielou evenness index (J) (Pielou, 1966), and the Margalef richness index (F) (Pielou, 1966). The effects of stand spatial structure on sapling diversity were analyzed using GAMs, following the methodology described above. The indices were calculated as follows:

Shannon-Wiener index:

Simpson index:

Pielou index:

Margalef index:

where S is the number of species, P_i is the proportion of individuals of species i relative to the total number of individuals, and N is the total number of individuals of all species.

The distributions of spatial structure parameters in 2011 and 2021 exhibited consistent patterns. For the uniform angle index, relative frequency increased initially and then decreased with higher levels, approximating a normal distribution. Across the sample plots, values of the uniform angle index generally fell within the range of random distribution, with more than half of the trees randomly dispersed throughout the study area. In contrast, very uniformly distributed trees were the least frequent. For the mingling index, low mingling levels consistently occurred at low frequencies, whereas relative frequency increased with higher mingling levels. In both years, the stand exhibited a strongly mixed state, indicating high species heterogeneity. The frequency distribution of the five dominance levels was relatively balanced, with overall dominance reflecting an intermediate state. This pattern suggests a neutral competitive balance among trees and stable equilibrium in individual competition dynamics (Table 3; Figure 1).

The spatial structure distribution of the stand exhibited a high degree of similarity between 2011 and 2021. In both years, the relative frequency within each dominance class consistently increased with higher levels of mingling, reaching its peak under complete mingling conditions. The maximum frequency occurred when the dominance was 0 and the mingling degree reached 1.00. Under zero mingling, the relative frequencies of trees across different competitive states were nearly identical, whereas higher mingling levels revealed more pronounced variations in frequency distribution among dominance classes (Figure 2).

For the uniform angle index, frequencies generally increased with greater mingling intensity, except under conditions of absolute uniformity. Within each mingling level, uniform angle index frequencies followed a normal distribution. The highest frequencies occurred under strong mingling combined with random spatial distribution. In contrast, frequencies were close to zero when absolute uniformity was paired with zero mingling, indicating that non-mingled trees were exceedingly rare in highly uniform stands (Figure 2).

Notably, when the uniform angle index was held constant, relative frequencies showed only minor variation across dominance classes. In contrast, within the same dominance level, uniform angle index frequencies exhibited a symmetrical decline from a central peak at random distribution, decreasing toward both clumped and uniform extremes. Randomly distributed trees constituted the majority of the stand, whereas absolutely uniform distribution was the rarest spatial pattern, underscoring its structural scarcity in the forest (Figure 2).

In both 2011 and 2021, the trivariate distribution of spatial structure showed that most trees within the same mingling and dominance classes exhibited random spatial patterns. Under identical mingling and uniform angle index conditions, relative frequency did not display a consistent trend with changes in dominance. However, at the same dominance and uniform angle index, frequencies increased gradually with rising mingling intensity. In both years, the highest relative frequencies occurred under random distribution combined with intermediate dominance and very high mingling intensity (Figure 3).

Elevation emerged as a key factor significantly influencing three spatial structure parameters. Overall, dominance increased with elevation, reaching its highest values at mid-slopes and hilltops. This pattern suggests that species size differentiation gradually decreased with elevation: trees in depressions exhibited greater dominance, whereas size differentiation was weaker at mid-slope and summit positions. A plausible explanation is that intense competition for water and light in rocky, high-elevation areas promotes the concentration of large-diameter individuals in depressions, while smaller trees persist at higher elevations, thereby reducing overall dominance. The pronounced variation in dominance across elevation gradients confirms elevation as a primary driver of size differentiation.

Elevation significantly influences plant distribution and growth both directly and indirectly by regulating key a factors such as light, temperature, water availability, and soil conditions (Zhao et al., 2005). The uniform angle index declined gradually with elevation before rising again, reaching its lowest values at 200–250 m. This pattern indicates more dispersed tree distributions in depressions and mid-slopes, but increasingly clustered patterns at hilltops. Elevational shifts in species richness are accompanied by intensified interspecific competition, self-thinning, and reduced intraspecific clustering. At higher elevations, increased rock exposure amplifies competition for light and heat, forcing trees into clustered distributions under harsh habitat conditions. Conversely, lower-elevation depressions and mid-slopes offer more favorable growth conditions, promoting dispersed tree distributions.

Elevation also strongly influences species segregation. The mingling index increased continuously in depressions, remained relatively stable at mid-slopes, and rose again at hilltops. Spatial segregation was lowest in depressions, moderate at mid-slopes, and highest at hilltop and bealocks. Areas with high segregation were characterized by strong species heterogeneity, elevated tree diversity, and more stable forest structures, often accompanied by a higher proportion of dominant species. This pattern likely reflects optimal light and hydrothermal conditions at higher elevations, which maximize resource availability, coupled with intensified intraspecific competition and self-thinning. However, Lv et al. (2023) reported that mingling increased with elevation while the dominance index and uniform angle index were unaffected, a discrepancy potentially attributable to karst habitat heterogeneity, where complex terrain alters structural parameter responses to elevation.

Slope and aspect significantly influence solar radiation, temperature, and precipitation distribution, thereby shaping the spatial arrangement and distribution of trees. Variations in slope gradient affect soil moisture content, which in turn regulates tree distribution. Gentle slopes exhibit superior soil water retention and nutrient preservation compared to steep slopes. As slope increased, the uniform angle index declined, reaching its minimum at approximately 40°, and then rose again beyond this threshold. Slope was identified as the primary driver of these changes in the uniform angle index. When the slope approached 40°, tree distribution tended tree distribution tended toward randomness, coinciding with optimal soil water-holding capacity and nutrient availability. In contrast, other slope ranges exhibited weaker soil conditions, resulting in more clustered tree distributions. Lin et al. (2021) reported that steep slopes are particularly susceptible to soil erosion due to shallow soil layers, reduced moisture, and lower nutrient availability, all of which constrain tree growth. In the northern tropical karst region, the steepest slopes are located at mid-hilltops, where prolonged solar exposure and extreme aridity result in severe soil moisture deficits, further limiting tree establishment.

Aspect primarily modulates light availability for trees, indirectly shaping heterogeneous microclimates in terms of temperature and humidity, which in turn affect tree distribution. With increasing aspect values, mingling declined, likely due to differences in sunlight exposure duration and humidity across aspects. In the northern tropical karst seasonal rainforest, hilltops receive direct solar radiation throughout the day, particularly in summer, when temperatures can exceed 60°C. Combined with high rock-bareness rate and intense evaporation, these conditions create extremely arid environments where only drought-tolerant plant species can persist. As a result, species segregation is markedly reduced on hilltops, with neighboring trees often belonging to the same species, forming simplified community structures. In contrast, enclosed depressions such as dolines buffer strong winds and accumulate rainfall, generating relatively humid microclimates. Under these conditions, species segregation is higher than on hilltops, supporting greater structural and compositional complexity.

Rock-bareness generates diverse microenvironments, creating high environmental heterogeneity while exerting strong filtering effects that significantly influence plant growth and distribution. In the 15-ha forest plot, the rock-bareness rate was approximately 95% near hilltops, 85% on mid-slopes, and only 10% at the valley bottom. As rock-bareness increased, the uniform angle index exhibited an upward trend, likely because limited soil depth and poor water retention in high rock-bareness areas forced species to cluster to access nutrients and moisture. In contrast, mingling declined with increasing rock-bareness, reflecting reduced species diversity under harsher conditions.

The hilltop extreme habitat, characterized by harsh, arid conditions, filtered out less-adapted species, leaving only drought-tolerant specialists such as Sinosideroxylon pedunculatum, Diospyros siderophylla, Pistacia weinmannifolia, and Boniodendron minius which are primarily distributed on mountain tops. In contrast, species such as Erythrina stricta, Vitex kwangsiensis, Sterculia monosperma, and Ficus hispida are mainly restricted to depressions with low rock-bareness rates. Consequently, high rock-bareness at the hilltops reduced species segregation, mitigating adverse effects and enhancing survival among drought-tolerant trees. Li et al. (2022) reported that in karst regions of southwestern China, rocks altered tree size but not species diversity, suggesting a neutral role in maintaining diversity. This finding is inconsistent with our results, which may reflect the presence of rare species in high-elevation barren zones, where complex microhabitats can support higher diversity. Conversely, in highly barren areas where nutrient scarcity constrains tree establishment, diversity tends to decline. Overall, variation in rock-bareness drives differences in plant light-resource utilization and ultimately shapes tree growth and distribution patterns (Chen et al., 2014).

The increase in mingling with rising mean diameter at breast height is consistent with previous studies (Pommerening and Uria-Diez, 2017; Wang et al., 2021, 2018). The mingling-size hypothesis, a spatial extension based on the Janzen-Connell effect and the population protection effect, posits that large trees are often surrounded by smaller, heterospecific neighbors and that species mixing is positively correlated with size inequality. Our findings support this hypothesis, suggesting its applicability to karst natural forests. Generally, large trees are associated with lower diversity but higher mingling compared to small trees. However, disturbances, management practices, and other ecological mechanisms may reduce mingling as tree size increases. Density-dependent processes can further modulate neighborhood species segregation at the community level. For small trees, dispersal limitation and environmental filtering often promote conspecific clustering, whereas large trees tend to display higher mingling.

To sustain growth, large trees compete intensely for resources, leading to self-thinning and long-term reductions in the number of conspecific neighbors. The resulting mortality creates microsites favorable for heterospecific colonization, thereby promoting mingling (Li et al., 2020). However, in our study, no significant difference in mingling was observed between large- and small-diameter trees, possibly reflecting the influence of stand developmental stage or local environmental conditions (Yang et al., 2019). Contrasting evidence from Knysna forests shows a negative correlation between mingling and size inequality, indicating that high spatial mingling need not associate with large trees (Pommerening et al., 2024). The positive mingling and diameter at breast height observed here may therefore reflect context-dependent dynamics, although such trend are reported in most stands.

Mingling decreased with increasing tree abundance, consistent with observations in the Nanpanjiang Basin (Pommerening et al., 2024). Rare species are more likely to be surrounded by heterospecific neighbors, whereas abundant species exhibit lower mingling, likely due to dispersal limitation. Species with short dispersal distances tend to form clumped distributions, increasing conspecific encounters and reducing mingling, a pattern supported by previous studies on dispersal constraints (Ma et al., 2024). As a result, greater species diversity enhances tree heterogeneity and promotes more uniform distributions.