Spatial autocorrelation shapes liana distribution better than topography and host tree properties in a subtropical evergreen broadleaved forest in SW China

Lianas are an important component of subtropical forests, but the mechanisms underlying their spatial distribution patterns have received relatively little attention. Here, we selected 12 most abundant liana species, constituting up to 96.9% of the total liana stems, in a 20‐ha plot in a subtropical evergreen broadleaved forest at 2472–2628 m elevation in SW China. Combining data on topography (convexity, slope, aspect, and elevation) and host trees (density and size) of the plot, we addressed how liana distribution is shaped by host tree properties, topography and spatial autocorrelation by using principal coordinates of neighbor matrices (PCNM) analysis. We found that lianas had an aggregated distribution based on the Ripley's K function. At the community level, PCNM analysis showed that spatial autocorrelation explained 43% variance in liana spatial distribution. Host trees and topography explained 4% and 18% of the variance, but less than 1% variance after taking spatial autocorrelation into consideration. A similar trend was found at the species level. These results indicate that spatial autocorrelation might be the most important factor shaping liana spatial distribution in subtropical forest at high elevation.


| INTRODUC TI ON
Knowledge on the ecological requirements determining plant distributions is a precondition for profitable and sustainable forest management and for forest conservation under current and future environmental conditions (Walthert & Meier, 2017). In contemporary forest ecology, the study of trees has remained the central focus, studies assessing the responses of liana species distribution to the environment are rare (but see Liu et al., 2021). Lianas are an important component of forest ecosystems (Schnitzer & Bongers, 2002). They compete with host trees for light in the canopy and water and nutrient resources belowground (e.g., Dillenburg et al., 1993;Schnitzer et al., 2005). As such, increases in lianas may alter forest structure and composition, decreasing tree growth, survival, and reproduction, resulting in reduced carbon sequestration and ecosystem functioning (e.g., Schnitzer et al., 2000;van der Heijden et al., 2015;Venegas-González et al., 2020). However, early research on lianas is largely focused on tropical forests (e.g. Schnitzer, 2018;Schnitzer & Bongers, 2002;van der Heijden et al., 2015). Subtropical forests have higher average net ecosystem productivity compared to tropical and temperate forests, accounting for 8% of the global forest net ecosystem productivity (Yu et al., 2014). Cascading effects from increasing lianas would potentially decrease subtropics' productivity; hence, the important role of lianas in subtropics is now gaining more attention (e.g. Ceballos et al., 2020;Yang et al., 2018;Yuan et al., 2009). Here, we evaluate the distribution patterns of the liana community in a subtropical forest in China.
Studies of subtropical lianas have made strides in identifying their composition and diversity, but little is known about the mechanisms that drive patterns of distribution. Liana distribution can be affected by abiotic and biotic factors, such as topographic features and forest structures (e.g., Kusumoto et al., 2013;van der Heijden & Phillips, 2008). For instance, the spatial distribution of lianas was controlled mainly by topographic gradients in a Japanese evergreen broadleaved forest, where most liana species had preferences for concave topographies (i.e., valley habitats, Kusumoto et al., 2013).
Other abiotic factors, such as light conditions mediated by canopy openness and tree-fall gaps (Malizia et al., 2010;Yuan et al., 2009), soil water, nutrients, and pH Malizia et al., 2010;Yang et al., 2018), and disturbance (Rice et al., 2004) also play important roles in liana abundance and distribution in subtropical forests. Lianas rely on host trees to ascend to the canopy with the help of stem twiner, tendril, root, and hook climbers (Putz, 1984). Early research (e.g., Malizia et al., 2010;Yang et al., 2018;Yuan et al., 2009) suggests that there might be a close relationship between liana abundance and host tree properties such as tree density, diameter, and height. For example, larger lianas tended to climb larger host trees in a subtropical forest in Taiwan . However, some studies have shown that host trees impose little or no influence on liana abundance and distribution in subtropical forests (Kusumoto et al., 2013;Malizia et al., 2010;Yang et al., 2018, but see Yuan et al., 2009. Liana species' spatial distribution can also result from dispersal limitation and demographic stochasticity (e.g., Aiba et al., 2012;Clark et al., 2018;Hubbell, 2001;Legendre et al., 2009). Early research has shown that most species, especially lianas, are spatially aggregated due to their short-ranged dispersal (Clark et al., 2018;Condit et al., 2000). Nevertheless, no studies have distinguished the relative contributions of environment and spatial autocorrelation to variability in liana distribution in the subtropical forest.
Research on subtropical lianas has been conducted in South America (e.g., Campanello et al., 2007;Ceballos et al., 2020;Malizia et al., 2010) and Asia (e.g., Kusumoto et al., 2008;Yang et al., 2018), with nearly all sites occurring at low elevations (below 1100 m) except for Ailaoshan (>2000 m, Yuan et al., 2009) in China. However, the distribution patterns of subtropical lianas are still unclear, which hampers the further understanding of the liana community across the world. In this study, our goal was to quantify the spatial distribution patterns of lianas in a 20-ha forest dynamic plot in Ailaoshan and how liana spatial distribution was related to biotic and abiotic factors. Specifically, we addressed how liana distribution is affected by host tree properties (density and size), topography (convexity, slope, aspect, and elevation), and spatial autocorrelation.

| Study site
This study was conducted in the Ailaoshan subtropical forest, a representative of the subtropical evergreen broadleaved forest in South China (Wu, 1980), possessing relatively mild winter and more rainy and foggy days compared to the eastern part of the subtropics in China (Zhang et al., 2013). We conducted the liana census in the Ailaoshan 20-ha forest dynamics plot (500 × 400 m), which was established based on the protocols of the Center for Tropical Forest Science in Ailaoshan National Nature Reserve, Yunnan Province, SW China (101.03°-101.03°E, 24.53°-24.54°N; 2472 to 2628 m a.s.l.; Wen et al., 2018). The plot was divided into 500 quadrats (20 × 20 m), and all free-standing woody plants with a diameter at breast height (DBH) ≥1 cm were mapped, tagged, measured, and identified to species level .
The mean annual precipitation is 1778 mm (86% in the rainy season from May to October) -2011Wen et al., 2018), and the mean relative humidity is 86%. The mean annual temperature is 11.3°C, with a mean temperature of 5.7°C in the coldest month (January) and 15.6°C in the warmest month (July). At a 0-20 cm depth, soil organic matter content was 170.4 ± 23.2 g kg −1 , total N, P, and K concentrations were 5.64 ± 0.77, 0.97 ± 0.09, and 9.46 ± 0.92 g kg −1 , available N, P, and K concentrations were 51.17 ± 7.70, 11.10 ± 2.78, and 185.39 ± 33.15 mg kg −1 , and pH was 4.24 ± 0.06 (data provided by the National Forest Ecosystem Research Station at Ailaoshan).

| Liana census
From October 2015 to January 2016, all rooted lianas with the diameter ≥1 cm were tagged, mapped, measured, and identified to species level following standard methods (Gerwing et al., 2006;Schnitzer et al., 2008). Briefly, we labeled and mapped the largest "principal stem" and additional multiple stems and then measured the diameter 130 cm from the rooting point at the soil surface. When lianas branched below 130 cm (but ≥40 cm from the roots), we measured the diameter 20 cm below the branching point. Because multiple stems seem to be functionally equivalent in their dynamics and effects on trees, we treated them as separate individuals in the subsequent statistical analyses (Schnitzer and Carson, 2001). We classified climbing mechanisms into four categories: stem twiner, tendril, hook, and root climbers (Putz, 1984). In total, 3725 liana stems, including multiple stems, were recorded in the 20-ha plot, representing 24 species and 15 families (Table 1). We selected the 12 most abundant species for the study (abundance >1 stem ha −1 ), constituting up to 96.9% (3609 stems) of the total liana stems and 97.5% (8.83 m 2 ) of the total basal area. We quantified the numbers of stems (including multiple stems) per species in each of the 500 quadrats.

| Topographic variables
For each of the 500 20 × 20 m quadrat, we calculated elevation as the mean elevation of the four corners of the quadrat. We calculated slope as the average angular deviation of each of the four triangular planes formed by connecting three corners at a time with respect to the horizontal line. Convexity was the difference between the average elevation of the focused quadrat and the average elevation of the eight neighbor quadrats (Song et al., 2018).
For the edge quadrat, convexity was calculated as the difference between the elevation of the center of the focal quadrat and the average elevation of the four corners. Aspect was calculated as follows 180 − tan −1 fy fx * 180 + 90 fy fx , where fx and fy are the elevation difference from east to west, north to south in the 20 × 20 m quadrat, respectively (Song et al., 2018).

| Statistical analyses
To detect clustering or dispersion trends, we used Ripley's K function to test the spatial distribution pattern of the 12 most abundant liana species (Ripley, 1976). Ripley's K function used with 999 Monte Carlo simulations and 95% confidence interval to analyze and compare the distribution patterns (regular, clustered, or random) of plant individuals. One of the most widely used methods to quantify the relative contribution of topography and spatial structure to species distributions is variation partitioning using topographic variables and Principal Coordinates of Neighbor Matrices (PCNM). PCNM is a type of Moran's eigenvector map and is closely related to Moran's index of spatial correlation (Dray et al., 2006). PCNM variables are used as spatial predictors in multivariate regression. In this analysis, pure effects of spatial structure are considered as dispersal limitation and other factors that are spatially structured but not related to topography. Most studies that distinguished the relative contributions of environment and spatial structure to species distribution have been conducted at the community level (Borcard & Legendre, 2002;Dray et al., 2006;Legendre et al., 2009). However, because plant species differ in their environmental requirements and dispersal limitation, each species may have unique responses to environmental and spatial variables.
We generated spatial variables using PCNM to characterize spatial structure at multiple scales (Borcard & Legendre, 2002;Dray et al., 2006;Legendre et al., 2009). The PCNM variables represent a spectral decomposition of the spatial relationships among sampled sites. First, we generated a matrix of Euclidean distance among 500 20 × 20 m quadrats. This matrix was then truncated using a defined threshold value. If the Euclidean distances were under the threshold, they were kept as measured, otherwise they were considered to be "large," and then this value will be set to four times the threshold value (Borcard & Legendre, 2002). We calculated the principal coordinate of this matrix and obtained 166 PCNM variables with positive Moran's I values. We then used these PCNM variables as spatial explanatory variables in the subsequent analyses.
We performed variation partitioning (Borcard et al., 1992) using a redundancy analysis to quantify the contribution of host trees (tree density and DBH), topography (elevation, slope, aspect, and convexity), and spatial variables to liana spatial distribution. To evaluate the species differences in liana spatial distribution, we also performed variation partitioning for each species based on multiple regression (Aiba et al., 2012). We separately performed forward selection for host trees, topography, and spatial variables and used selected variables for the final model to avoid overestimation of R 2 (Blanchet et al., 2008). This step was only performed at the community level and the same suite of variables was used for species-level regressions to minimize the risk of selecting superfluous variables through repetition of variable selection (Aiba et al., 2012). We then performed seven multiple regressions for each species to obtain the percentage of variation explained by one to three variables from the host trees, topography, and spatial variables. All statistical analyses were performed in the R programming environment v. 3.6.0 (R Core Team, 2019).

| RE SULTS
In total, 3725 liana stems including multiple stems (186 lianas ha −1 ) were recorded in the 20-ha plot, representing 24 species and 15 families ( Table 1). The differences in species distribution were related to differences in the topographic habitats in Ailaoshan sub-  semicordatatifolia species with total stems ≥20 in the 20-ha plot were spatially aggregated at most spatial scales (Figure 2).
Redundancy analysis showed that, at the community level, 4% of the variation in liana spatial distribution was explained by host tree properties (tree density and DBH), 18% was explained by topographic variables (convexity, slope, aspect, and elevation), and 43% was explained by spatial autocorrelation (Figure 3a). The effects of both host trees and topography were strongly spatially structured. After controlling for the effect of spatial autocorrelation, the contributions of the host trees and topography were 0% and 1%, F I G U R E 2 Plots of Ripley's K(r) versus scale r (m) for 12 liana species with more than 20 stems in the 20-ha forest dynamics plot. The solid line shows observed values, the dashed line denotes the complete spatial randomness, and the gray area represents the 95% confidence interval (Monte-Carlo = 999 times). Values above the gray area (totally random distribution curve) indicate a clumped distribution, within the gray area indicates a random distribution and below the gray area a regular distribution  (Figure 3a). Multiple regressions showed that partitioning of variance at the species level was similar to that at the community level (Figure 3b). Similarly, spatial autocorrelation explained the highest variation in spatial distribution at the species level. The contribution of spatial autocorrelation ranged from 5.0% to 22.0%, followed by topography (from 1% to 5.5%) and host trees (up to 2.4%).

| DISCUSS ION
In this study, we investigated liana spatial distribution patterns and their drivers in a subtropical evergreen broadleaved forest in SW China. We found that liana distribution in the Ailaoshan forest was spatially aggregated. Spatial autocorrelation explained most of the variation in liana distribution, with the number of host trees and topographic variables having little effect. This is consistent with other local-scale studies in subtropical forest (Malizia et al., 2010;Yang et al., 2018). Kusumoto et al. (2013) found that host tree availability (species richness and density) had little effect on liana species distribution. This finding does not rule out the importance of host trees but may be ascribed to canopy structure of the subtropical forest.
For example, the host tree size has strong influence on liana abundance in tropical forests (e.g., Homeier et al., 2010;Leicht-Young et al., 2010;Poulsen et al., 2017). In addition, tree canopy might indirectly affect the liana community by modifying light conditions in forests (van der Heijden & Phillips, 2008), and liana abundance and distribution are sometimes strongly correlated with forest canopy openness and height (Yuan et al., 2009). These features may facilitate liana increase, making them more abundant in places with better light conditions such as tree-fall gaps (Malizia et al., 2010;Yuan et al., 2009).
Topography had a relatively stronger impact than host trees on liana spatial distribution in our study. At the species level, R. longicuspis was more abundant at low elevations and flat habitat and K. heteroclita was more abundant in valleys. Our results are in line with reports from other subtropical forests which revealed that lianas have higher abundance and diversity in valley sites (Kusumoto et al., 2013;Yang et al., 2018) and low-elevation habitats . Even in tropical forests, lianas have high abundance and diversity in valleys (Putz & Chai, 1987;Romero-Saltos, 2011) and low-elevation sites (Addo-Fordjour & Rahmad, 2015). Valleys normally have high soil moisture and nutrients, which allow lianas to meet their high nutritional and moisture demands for rapid growth (Kusumoto et al., 2008;Malizia et al., 2010;Yang et al., 2018;Yuan et al., 2009).
Liana distribution in the Ailaoshan subtropical forest was strongly driven by spatial autocorrelation, with host trees and topography explaining very little of the variation in liana distribution after taking spatial autocorrelation into consideration at the community and the species levels. Similarly, spatial processes strongly influenced liana abundance and distribution in Barro Colorado Island, Panama, where the average dispersal distance of lianas is ca. 3 m, resulted in their aggregated distribution (Clark et al., 2018). In our study, 10 of 12 liana species are dispersed via zoochory, which tend to be more aggregated due to short-distance dispersal (Clark et al., 2018). Therefore, it is likely that dispersal limitation also shapes liana spatial patterns in subtropical forests, which deserves further investigation in the future.
A high proportion of the variation in liana abundance (56%) remained unexplained, mainly due to two reasons. First, we did not measure the nutrient status and soil moisture in the study plot, which could strongly affect liana spatial distribution (Kusumoto et al., 2013;Malizia et al., 2010;Yang et al., 2018;Yuan et al., 2009). For instance, soil P concentration was strongly and positively related to liana abundance in a subtropical montane forest in Northwest Argentina (Malizia et al., 2010). Furthermore, soil P, pH, and moisture have been found to be positively related to liana abundance in tropical forests in Malaysia (e.g., Addo-Fordjour et al., 2014), Ghana (Addo-Fordjour & Rahmad, 2015) and southwestern China . Second, dispersal limitation could be a major mechanism underlying liana species distribution (Clark et al., 2018).
However, dispersal traits such as seed mass were not included in the F I G U R E 3 Variation partitioning at the community (a) and species (b) levels for spatial abundance distribution for each of the 12 species with more than 20 stems in the 20-ha forest dynamics plot. See Table 1

| CON CLUS ION
This is the first study to distinguish the relative contributions of host trees, topography, and spatial autocorrelation on liana distributions in a subtropical evergreen broadleaved forest at high elevation in SW China. Our results showed that the spatial distribution of lianas in study area is strongly aggregated, a result of spatial autocorrelation rather than by host tree size and abundance or topographic variables. Therefore, to conserve plant diversity in subtropical forest ecosystems at high elevation, further investigation of liana diversity maintenance and community assembly in subtropical forest ecosystems is necessary.

ACK N OWLED G M ENTS
This study was financially supported by the National Natural

AUTH O R CO NTR I B UTI O N S
XLB and JLZ designed the study. XLB and HDW performed the experiments. XLB and QL analyzed data. XLB, QL, JLZ, DM, and SKG led the manuscript writing. All authors contributed intellectually to the manuscript, revised it critically, and gave final approval of the version to be submitted.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are openly available in the Dryad Digital Repository: https://doi.org/10.5061/dryad. tqjq2 bw10 (Bai et al., 2021).