The Interaction between Climate and Soil Properties Influences Tree Species Richness in Tropical and Subtropical Forests of Southern China

: The climate and soil properties are major determinants of plant growth and forest community assembly across diverse biomes. However, the contribution of climate and soil on species diversity in tropical and subtropical forests remains controversial. Therefore, this study aims to evaluate the effects of soil and climate on tree species richness using survey data across 495 tropical and subtropical forest plots in Southern China. The selected predictors were categorized as local plot characteristics, climate, and soil factors, and their relationship with tree species richness was modeled using negative binomial generalized linear models. The results revealed that the considering of the interaction between climate and soil properties considerably improved the goodness − of − fit of these models. The individual effects of climate and soil factors had weak relationships with species richness, accounting for 3.61% and 5.77% of the overall 58.9% explained variance in species richness, respectively. Instead, the interaction between climate and soil properties explained most of the statistical variation in tree species richness (84.34% of the overall 58.9% explained variance). The results highlight the importance of soil and climate interactions on tree diversity in tropical and subtropical mature natural forests and their incorporation into biodiversity assessment models to enhance the prediction of community change and responses to climate change.


Introduction
Environmental changes are altering the pattern and the assembly processes of biodiversity, leaving it unclear whether the persistence of biodiversity patterns would be impacted positively or negatively by these environmental changes [1][2][3][4].Therefore, clear understanding of abiotic factors underlying species diversity across communities is vital for monitoring community stability and structural evolution under changing climate, thus designing appropriate conservation and management measures [5,6].
Biomes and forest community distribution around the globe match variation in climate and soil properties.Notably, as the crucial factor that has a direct impact on the survival, growth, and distribution of plants, climate has always been the primary factor to explain the diversity patterns of plant species in primary forests [7][8][9].Studies undertaken at continental or global scales revealed the predominant effect of climatic variables in explaining

Study Area and Data Source
This study was undertaken in multiple vegetation types across Guangxi, China (Figure 1).The forest survey data of 495 plots in mature natural forests were collected from two major databases including "The vegetation of Guangxi [47]" and "Vegetation of Guangxi (2 volumes) [48]".These survey data included the geographic name, latitude, longitude, elevation, area, species composition, and abundance (or relative abundance) of each plant species within each plot.The latitude and longitude of 13 plots (2.62%) without precise geographic coordinates were approximated from the locality name (village) using Google Earth [49].The area of the plots ranged from 100 to 1500 m 2 with 96.57% being smaller than 1000 m 2 , and of which 47.47% were 400 m 2 and 27.47% were 600 m 2 , respectively.
Forests 2024, 15, x FOR PEER REVIEW 3 of 16 mainly explains the variation in the tree species richness in mature natural forests in this region.

Study Area and Data Source
This study was undertaken in multiple vegetation types across Guangxi, China (Figure 1).The forest survey data of 495 plots in mature natural forests were collected from two major databases including "The vegetation of Guangxi [47]" and "Vegetation of Guangxi (2 volumes) [48]".These survey data included the geographic name, latitude, longitude, elevation, area, species composition, and abundance (or relative abundance) of each plant species within each plot.The latitude and longitude of 13 plots (2.62%) without precise geographic coordinates were approximated from the locality name (village) using Google Earth [49].The area of the plots ranged from 100 to 1500 m 2 with 96.57% being smaller than 1000 m 2 , and of which 47.47% were 400 m 2 and 27.47% were 600 m 2 , respectively.

Environmental Predictors
Environmental predictors used were divided into three categories including local plot characteristics, climate variables, and soil variables.The local plot characteristics included the area of the plot, longitude, latitude, elevation, forest type, and habitat type.Due to the power function relationship between species richness and area, the area of the plot was transformed into a logarithmic form for facilitating the regression [50,51].The elevation data were extracted from the SRTM (Shuttle Radar Topography Mission) elevation map that has been downloaded from the Worldclim website (https://www.worldclim.org/(accessed on 9 September 2019)), at a 30 arc−second resolution (about 1 km at the equator).The forest types of the plots were categorized into seven forest types based on the World Widelife Fund for Nature (WWF) forest classification standards and the types reported by the survey operator, namely tropical forest (TF), evergreen broad−leaved forest (EBF), deciduous and evergreen broad−leaved mixed forest

Environmental Predictors
Environmental predictors used were divided into three categories including local plot characteristics, climate variables, and soil variables.The local plot characteristics included the area of the plot, longitude, latitude, elevation, forest type, and habitat type.Due to the power function relationship between species richness and area, the area of the plot was transformed into a logarithmic form for facilitating the regression [50,51].The elevation data were extracted from the SRTM (Shuttle Radar Topography Mission) elevation map that has been downloaded from the Worldclim website (https://www.worldclim.org/(accessed on 9 September 2019)), at a 30 arc−second resolution (about 1 km at the equator).The forest types of the plots were categorized into seven forest types based on the World Widelife Fund for Nature (WWF) forest classification standards and the types reported by the survey operator, namely tropical forest (TF), evergreen broad−leaved forest (EBF), deciduous and evergreen broad−leaved mixed forest (DEBMF), deciduous broad−leaved forest (DBF), coniferous and broad−leaved mixed forest (CBMF), and evergreen coniferous forest (ECF).Additionally, the habitat type of each plot was identified as karst or non−karst landscapes based on the karst spatial distribution map obtained from the National Earth System Science Data Center (https://www.geodata.cn/main/(accessed on 9 April 2018).
The climatic data were obtained from the China ground−based average dataset of daily meteorological records from 1981 to 2010, recorded by 2160 base, standard, and ordinary ground−based meteorological observatories, which can be downloaded from the China Meteorological Data Service Center (https://data.cma.cn/(accessed on 8 March 2019)).The dataset included six daily variables: average temperature, average maximum temperature, average minimum temperature, average vapor pressure, average precipitation, and average wind speed [53].
To generate more biologically meaningful variables, we calculated 19 bioclimatic variables and the potential evapotranspiration (PET) for each plot as the climate variables.First, we extracted the meteorological data of the 75 observatories distributed in Guangxi.Then, 19 bioclimatic variables and PET of observatories were calculated according to the descriptions of bioclimatic variables on the website of Global Climate Data (https: //worldclim.org/(accessed on 9 September 2019)) and the Penman-Monteith equation recommended by the Food and Agriculture Organization of the United Nations (FAO), respectively.After detrending these 19 bioclimatic variables and PET of observatories using a quadratic polynomial, we sequentially applied 8 kriging models to computerize the semi−variogram of residual distribution for these detrended variables.By visual inspection and parameter assessment, the spherical kriging model with the smallest sum of squared residuals is selected as the optimal one.Finally, we calculated the values of the 19 bioclimatic variables and PET at each location of 495 plots by this spherical regression-kriging approach using the "mgcv" and "gstat" packages in R [53] (see Table S1 in Supplementary Materials for more details).These variables are widely used in species distribution modeling and related ecological modeling [39,41,51,54,55].
Many intercorrelations were observed between the 20 climate variables and 27 soil variables, respectively.Therefore, we excluded the variables with a correlation coefficient greater than 0.6 using a canonical correlation test.Finally, 6 of 20 climate variables (Bio3, Bio7, Bio8, Bio12, Bio15, PET) and 15 of 27 soil variables (AN, BD, CA, CEC, CL, DH, H, K, MG, POR, SW1, TK, TN, TP, WH) were retained.All variables used in this study were standardized by the z−scored method (mean−centered and divided by the standard deviation).

Modeling Procedure
The negative binomial generalized linear model with a log link function has been used in this study because the species richness values of all plots showed a negative binomial distribution (see Figure S1 in Supplementary Materials).The species richness is the dependent variable, and local plot characteristics, climate, and soil variables are predictors.Our interest here is to explore the individual effects of local plot characteristics, climate, soil, and the joint effects of climate and soil variables on shaping the species diversity pattern of trees in tropical and subtropical forests.The quality of models was based on five combinations of predictors including (1) only the local plot characteristics; (2) local plot characteristics and climate variables; (3) local plot characteristics and soil variables; (4) local plot characteristics, climate, and soil variables; (5) local plot characteristics, climate, soil variables, and the interactions between climate and soil variables.After constructing the full models of all combinations of predictors, the stepwise model selection method was used to simplify the model based on the Akaike Information Criterion (AIC).
To avoid the overfitting, only one plot was randomly selected from each 30 arc−second grid, so that 293 of all 495 plots were left for modeling.The relative independence of the distribution of plots could reduce the bias caused by the spatial autocorrelation of an aggregated distribution of plots [56][57][58].Plot selection and model fitting were repeated 999 times, and then the average and standard deviation of the model coefficients were calculated based on the 999 repeats.The relative importance of the predictors was calculated based on the percentage of variance explained by variables explained using the ratio between the absolute values of their standardized regression coefficients and the sum of all standardized regression coefficients from the predictors [59].This method is similar to a variance partitioning analysis [59][60][61].
The data processing was performed in ArcGIS (version 10.2) and R software (version 4.3.2;https://www.r-project.org/(accessed on 11 December 2023)).The negative binomial generalized linear model was constructed using the function "glm.nb" in the "MASS" R package.

Performance of Five Combinations of Predictors
The prediction of the best selected model for each combination of predictors is shown in Figure 2. The frequency distributions (left panel in Figure 2) and distribution maps (right panel in Figure 2) of observed and predicted species richness of trees among 293 plots both showed that the individual effects of local plot characteristics, climate, and soil factors had weak relationships with species richness.The local plot characteristics only explained 22.6% of the total variation in species richness among plots.The independent consideration of climate (Figure 2C,D) and soil variables (Figure 2E,F) increased the model power (R 2 ) by 3.7% and 9.2%, respectively.The addition of both climate and soil variables into the model (Figure 2G,H) increased the model power by 11.1%.The species richness predicted from the above models was significantly different from the observed value.However, the consideration of the climate, soil variables, and their interactions in the model significantly improved the power of the model with a total R 2 of 58.9%, which was 36.3% higher than the model using the local plot characteristics only, thus representing the best fitting model among all the combinations of predictors (Figure 2I,J).Furthermore, when considering the interaction between climate and soil properties, the predicted species richness (the radius of the circle in the right panel in Figure 2) exhibited a wider range of variation.Although the joint effects of individual climate and soil variables may be complex, the overall prediction results are consistent.
The specific contribution of each variable to the variation in local species richness could be quantified by comparing the difference in the ratio between the absolute values of their standardized regression coefficients and the sum of all standardized regression coefficients from the predictors.Based on the best fit model, the individual effects of local plot characteristics, climate, and soil variables were 6.28%, 3.61%, and 5.77%, respectively (Figure 3).However, the interactions between climate and soil variables explained 84.34% of the overall 58.9% explained variance in species richness among plots.plot characteristics, climate, and soil variables were 6.28%, 3.61%, and 5.77%, respectively (Figure 3).However, the interactions between climate and soil variables explained 84.34% of the overall 58.9% explained variance in species richness among plots.

Effects of Predictors
The trend of species richness was positively correlated with the longitudinal gradient (  β = 0.52, p = 7.4 × 10 −4 ), but was not significant with the latitudinal gradient (Figure 4 and Table S2 in

Effects of Predictors
The trend of species richness was positively correlated with the longitudinal gradient ( β = 0.52, p = 7.4 × 10 −4 ), but was not significant with the latitudinal gradient (Figure 4 and Table S2 in the Supplementary Materials).The plot area ( β = 0.21, p = 9.1 × 10 −13 ) and habitat type ( β = −0.55,p = 0.05) significantly influenced the species richness across different plots.The effect of specific climate variables on local species richness may differ.The temperature annual range (Bio7, β = −0.41,p = 0.02) and the mean temperature of the wettest quarter (Bio8, β = −0.38,p = 0.02) were negatively associated with local species richness.Moreover, CEC ( β = 0.49, p = 1.1 × 10 −4 ), H ( β = −0.23,p = 0.04), MG ( β = −0.27,p = 4.7 × 10 −3 ), and TK ( β = 0.26, p = 0.05) within soil variables also significantly affected the regional variation of local species richness of trees.The potential evapotranspiration (PET) did not emerge as a significant independent variable in our study although this variable was associated with species richness in some research.The individual effects of climate and soil variables were weak, but almost (98.89%)all interactions between climate and soil variables significantly affected the local species richness.= −0.27,p = 4.7 × 10 −3 ), and TK (  β = 0.26, p = 0.05) within soil variables also significantly affected the regional variation of local species richness of trees.The potential evapotranspiration (PET) did not emerge as a significant independent variable in our study although this variable was associated with species richness in some research.The individual effects of climate and soil variables were weak, but almost (98.89%)all interactions between climate and soil variables significantly affected the local species richness.

Discussion
Our research suggests that a complex combination of climate, soil, and geographic variables significantly affected local species richness in the present region with high habitat heterogeneity.In this study, the individual effects of soil and climate variables were unable to explain the variation in species richness of trees.However, the interaction of soil and climate variables had a significant impact on the variation of species richness, which made the explanation rate increase by 25.2% relative to the independent effects of predictors.This is in line with previous studies that indicate that the diversity pattern is affected by both climate and soil factors, especially the interaction between climate and soil, which may dominate the composition and structural characteristics of local forests [10,29,[62][63][64].There is ample evidence to suggest that environmental factors such as climate and soil play a significant role in driving the patterns of species geographical distribution, species and functional diversity, and large−scale variation in plant traits across spatial scales, in which the joint effects of factors dominate the variation compared to individual effects [13,14,29,[65][66][67][68].These findings could be interpreted as evidence highlighting the importance of both joint and individual effects of climate and soil variables in shaping and maintaining the plant species diversity patterns across regional and global scales.
In recent years, the effects of climate and soil factors on diversity patterns of plants have received increasing attention [66,69].The effects of the interactions of climate and soil on species diversity patterns remains a complex puzzle.It is well known that climate is an important driver of forest structure and composition and can affect the soil property as well [29,70,71].Climate determines the amount of water and temperature associated with the soil weathering and the nutrient release.Under warm and humid climate conditions, soil microbial activity is faster than those under cool and dry climate conditions, stimulating the soil organic matter decomposition, thus facilitating the coexistence and growth of more plant species [72][73][74].Therefore, the interdependence between climate and soil can significantly alter species composition and diversity in natural ecosystems [75,76].For example, Palpurina et al. showed that the relationship between the species richness and soil pH value in dry steppes in eight regions of northern Eurasia was indirectly influenced by the precipitation [77].Ding and Eldridge suggested that climate and plants regulate the spatial variation in soil multifunctionality in eastern Australia [78].Moreover, the physical and chemical properties of soil also affect micro−climate at fine scales.Moist or dense soil is more likely to reduce the heat dissipation and thus stabilize the ambient temperature [79].Complex interdependent relationships exist between soil, climate, and species diversity.Advances in multivariate analyses such as structural equation modeling (SEM) may help with explaining these interactions.
In addition, one potential reason of the high relative effect of the interaction of climate and soil variables on species richness in primary forests may be the positive correlation between species diversity of plants and soil organic carbon (SOC) storage [80][81][82][83].The litter production and decomposition rate determine SOC storage and both are affected by air temperature and soil moisture [84,85].Although SOC−related variables were not used in this study, there are complex interactions between SOC and nitrogen and phosphorus, which are two abundant elements in soil [86,87].For example, soil nitrogen fixation can further affect the associated carbon-nitrogen cycling processes, including mycorrhizal symbiosis and litter decomposition [88].Under the high-phosphorus condition, the phosphatase released by soil microorganisms may be mainly used for the purpose of obtaining carbon sources [89].Soil is an important part of the carbon cycling processes.Differences in local soil physical and chemical properties can affect the emission of carbon−containing gasses by the soil at a small scale.This "climate-soil microenvironment" plays an important role in the growth and reproduction of local tree species [63].However, there is a two−way interaction between soil organic matter and plants.Higher plant diversity also directly increases the rhizosphere carbon input of the microbial community, thereby increasing microbial activity and carbon storage [90].The reason for selecting the data of the first soil layer in this study is the barrenness of the soil formed by the unique geological characteristics and hydrological structure in the karst areas.The physical and chemical properties of soil in karst areas are more variable and directly affected by bedrock types and complex topography because of the underdeveloped soil texture and shallower soil depths of karst areas [91,92].Karst habitats exhibit strong spatial heterogeneity in soil nutrients, soil properties, and water availability across the topographic position [93].Therefore, the large variation in soil nutrients can strengthen the correlation between vegetation and habitat conditions in karst areas.Forest survey and soil property measurement in Guangxi and other tropical karst forests revealed significant correlations between species diversity and soil nutrient contents [94][95][96].However, the low soil phosphorus content and high soil pH in karst areas led to a decrease in the bioavailability of phosphorus [97].This may be a potential reason why the individual effect of total soil phosphorus concentration on species richness is not significant in this study.
The species richness of forest trees also varies due to the unique historical and geographical characteristics of continental regions, independent of local climate and soil conditions.Qian et al. pointed out that species and phylogenetic diversity in East Asia are significantly greater than in eastern North America, and fewer species are contained in these forest plots in North America than that expected based on the statistical relationship between diversity and current local climate [98].The local species diversity of trees in primary tropical and subtropical forests is the result of large−scale and local−scale assembly processes.If the speciation and extinction would be affected by the specific historical and biogeographic events at a regional scale, then the traces of these specific events on the regional pattern of local species richness might be detected [51,99].Due to the continuous uplift of the Himalaya Mountains and the continuous occurrence of subsequent orogeny, the Hengduan Mountains and the Qinling Mountains in East Asia have formed an "Г" shape, forming the stable environment conditions conducive to the growth and reproduction of plants [98,100].We believe that the geo−historical processes potentially shape the geo−biological and diversity patterns of Guangxi.
The present study contains some limitations and we advise readers to take precautions when interpreting the results.For instance, most of the 495 plots used in this study were surveyed in the 1980s and 1990s, and a small number of plots were surveyed in the first decade of the 21st century.At that time, disturbance such as wildfire and over−logging may have affected the vegetation integrity and species composition within the dataset.Indeed, forest structure and species composition strongly differ over time and space [101], but we failed to consider the level of temporal variation in the analysis.Now, after decades of succession, the species composition of forest communities in different plots may have changed to varying degrees.Additionally, inequal plot size can result in some statistical issues.Indeed, the plot size varied from 100 to 1500 m 2 , of which nearly half of the plots were 400 m 2 , and more than 95% of the plots were less than 1000 m 2 .Rare species can hardly be inventoried using small plots, underestimating the species richness and weakening the correlation between the diversity pattern and environmental factors.
The potential evapotranspiration (PET) was retained in this study due to the significant correlation between local species richness and the interaction of PET and soil variables.However, the individual effect of PET did not show a significant correlation with local species richness.The role of variables such as PET in explaining the variation in diversity patterns varies across different studies.The specific evaluation of the role of certain variables at different scales is still necessary.For instance, PET was found to be more indicative of the effect of heat on local plant diversity than temperature in southern Africa [102], and it had a positive effect on explaining the phylogenetic beta diversity within 30 1 ha plots in Guangxi at a regional scale [53].However, it did not show a significant correlation when evaluating the impact of region effects on local species diversity within 47 forest plots worldwide [51].
Most environmental factors from a field survey always have complex, strong collinearity relationships, including the bioclimatic and soil factors [103,104].The generalized linear model that theoretically requires weak correlations among these factors is unrealistic.Therefore, the threshold of the Pearson correlation coefficient has become an important parameter that determines the number of associated predictor variables, which would be defined and excluded to adjust the number of used variables and their contribution.If fewer associated predictor variables are excluded because of the high threshold, the average contribution of a single variable would be relatively small, which in turn reduces the statistical significance of the model.This threshold is generally recommended to be 0.70 [105].In this study, after repeated evaluations, we found that the model with a threshold of 0.6 performed relatively well in both realistic and statistical terms.
Scale dependence plays a very important role in diversity research [40].Ecologists generally agree that the species diversity pattern at the small and medium scales is a balance between adaptation to environmental conditions and competition for resources [20,32,106,107].Many arguments on environmental filtering assume that the abiotic and biotic drivers of community assembly can be seen as independent effects, representing environmental filtering and interspecific/intraspecific interactions as the successive steps during the assembly processes [31,[108][109][110].However, in practice, these two factors interact concurrently and dynamically to drive community assembly [37,111].The intensity and direction of interspecific/intraspecific interactions could be strongly influenced by the abiotic conditions [112], even by the abiotic filtering at different life stages [113,114].Therefore, the effects of the interactions between climate and soil variables in this study probably include a general effect of abiotic processes and biotic processes that potentially affect by abiotic processes at the local scale.Consequently, this study advances our understanding of the patterns of plant species diversity at regional scales.

Conclusions
We modeled the relationship between local species richness and the predictors including the local plot characteristics, climate, and soil variables at a regional scale by the negative binomial generalized linear model.The results showed that the species richness of trees in natural mature forests was mainly affected by the interaction between climate and soil variables, but the individual effects of climate and soil variables were relatively low.We highlighted individual and joint effects of climate and soil on species richness of trees, which have been neglected so far.We claim that more accurate and informative census datasets are still needed to improve ecological analysis and modeling assessments, and to further understand the role of the abiotic filtering in shaping the diversity patterns of plant communities.

Supplementary Materials:
The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/f15081441/s1,Table S1: Meanings of the 19 bioclimatic variables; Table S2: The output of the negative binomial generalized linear model for the local species richness at a regional scale; Figure S1: The frequency distributions of observed species richness of trees among 293 plots for 999 random selections and the fitted negative binomial distributions (red curves).

Figure 1 .
Figure 1.The distribution map of the 495 plots in Guangxi Zhuang Autonomous Region, China.The color background represents the elevation distribution.

Figure 1 .
Figure 1.The distribution map of the 495 plots in Guangxi Zhuang Autonomous Region, China.The color background represents the elevation distribution.

Figure 2 .
Figure 2. The frequency distributions (left panel) and distribution maps (right panel) of observed and predicted species richness of trees among 293 plots by the negative binomial generalized linear model for 999 replicates under different combinations of predictors.The combinations of predictors in the figure are (1) only local plot characteristics (A,B); (2) local plot characteristics and climate variables (C,D); (3) local plot characteristics and soil variables (E,F); (4) local plot characteristics, climate, and soil variables (G,H); (5) local plot characteristics, climate, soil variables, and the interactions between climate and soil variables (I,J).In the left panel, the error bars on the red curves are

Figure 2 .
Figure 2. The frequency distributions (left panel) and distribution maps (right panel) of observed and predicted species richness of trees among 293 plots by the negative binomial generalized linear model for 999 replicates under different combinations of predictors.The combinations of predictors in the figure are (1) only local plot characteristics (A,B); (2) local plot characteristics and climate variables (C,D); (3) local plot characteristics and soil variables (E,F); (4) local plot characteristics, climate, and soil variables (G,H);(5) local plot characteristics, climate, soil variables, and the interactions between climate and soil variables (I,J).In the left panel, the error bars on the red curves are the standard deviation of predicted species richness among 999 replicates.In the right panel, for each point, the x coordinate of the center of the circle is the observed species richness, the y coordinate of the center of the circle is the average value of the predicted species richness of 999 replicates, the radius of the circle is the standard deviation of the predicted species richness of 999 replicates, and the legend is the abbreviation of the forest type.The values of R 2 are all adjusted R 2 .If the model predicted species richness perfectly, the points would fall on the diagonal lines (right panel).CBMF: coniferous and broad−leaved mixed forest; DBF: deciduous broad−leaved forest; DEBMF: deciduous and evergreen broad−leaved mixed forest; EBF: evergreen broad−leaved forest; ECF: evergreen coniferous forest; TF: tropical forest.

Figure 3 .
Figure 3. Importances of five combinations of predictors to explain the variance in species richness among plots for 999 replicates.All predictors were standardized to interpret parameter estimates on a comparable scale.

Figure 3 .
Figure 3. Importances of five combinations of predictors to explain the variance in species richness among plots for 999 replicates.All predictors were standardized to interpret parameter estimates on a comparable scale.

Figure 4 .
Figure 4. Average relative effects of local plot characteristics, climate, soil variables, and the interactions between climate variables and soil variables on species richness among plots for 999 replicates.The relative importance of each predictor is expressed as the percentage of the explained variance (left panel).The average parameter estimates (standardized regression coefficients) and the relevant average 95% confidence intervals are shown in the right panel.The right panel represents retained predictors of the best model selected by AIC.The adjusted R 2 of the averaged model and the p value for each predictor have been given as * p < 0.05; ** p < 0.01; and *** p < 0.001.(See TableS2in Supplementary Materials for detailed numerical results.)

Table S2 in
Supplementary Materials for detailed numerical results.)