Spatial patterns and interspecific associations among trees at different stand development stages in the natural secondary forests on the Loess Plateau, China

Abstract Quercus wutaishansea populations on the Loess Plateau are currently becoming more dominant in natural secondary forests, whereas Pinus tabulaeformis is declining. In the present paper, the diameter class (instead of age) was used to classify the different growth stages as juvenile, subadult, or adult, and the univariate function g(r) was used to analyze the dynamic changes in spatial patterns and interspecific associations in three 1‐ha tree permanent plots on the Loess Plateau, NW China. Our results suggested that the niche breadth changed with the development stage. The diameter distribution curve was consistent with the inverted “J” type, indicating that natural regeneration was common in all three plots. There was a close relationship between the spatial pattern and scale, which showed significant aggregation at small distances, and became more random as distance increased, but in the Pinus + Quercus mixed forests, the whole species were aggregated at distances up to 50 m. The degree of spatial clumping decreased from juvenile to subadult and from subadult to adult. The spatial pattern also differed at different growth stages, likely due to strong intraspecific competition. Associations among different growth stages were positively correlated at small scales. Our study is important to the understanding of the development of the Q. wutaishansea forests; thus, the spatial dynamic change features should be received greater attention when planning forest management and developing restoration strategies on the Loess Plateau.

be a consequence of many diverse factors, which include regeneration strategies such as seed dispersal, clonal propagation, soil conditions, and disturbances that affect the canopy structure (Larson et al., 2016;Noma, 1997;Yousef & Krzysztof, 2017). In contrast, regular spatial patterns are often attributed to local density-dependent influences, which include the mortality caused by species-specific natural enemies and competition from conspecific neighbors (Fibich, Leps, Novotny, Klimeš, & Těšitel, 2016;Tanner, Hughes, & Connell, 1994). The spatial patterns of trees may also be influenced by biotic and abiotic factors during stand development, and the appropriate disturbance-based silvicultural strategies could be designed to influence forest dynamic processes (Whitfeld et al., 2014).
Interspecific associations could indicate the spatial distribution relationship and functional dependency of species and their environment (Su et al., 2015), which is central to successional theory.
Complex and stable forest communities result from the replacement and establishment of different tree species during stand development, and changes in tree size may be informative in the study of developmental changes across forests because of competition among similarly sized individuals (Moeur, 1993;Salas, Lemay, Nunez, Pacheco, & Espinosa, 2006). Several models have implicitly assumed that interspecific associations are of primary importance in determining stand development patterns (Janík et al., 2016;Nakashizuka, 2001). The arrangement of tree species could provide significant information about forest dynamics processes and inter-and intraspecific interactions (Alekseev & Zherebtsov, 1995), including forest establishment, tree growth, species competition, plant reproduction, and mortality (Martíneza, Wieganda, Taboadab, Fernando, & Joséramón, 2010;Schleicher, Wiegand, & Ward, 2011).
The Loess Plateau is located in the northwestern China and is one of China's four large plateaus. This plateau is one of the birthplaces of ancient civilizations in China and features the most concentrated distribution and largest areas of loess on Earth. Most importantly, the Loess Plateau contains some of the largest areas of soil erosion and possesses one of the most vulnerable ecosystems in the world (Sun et al., 2015). Reforestation or natural forest regeneration can be an effective measure in controlling the soil erosion and thus in ameliorating the ecological environment. During reforestation or natural forest regeneration, internal competition among plants for limited resources has been proven to be the driving factor for forest succession. Quercus wutaishansea forests are the dominant or natural climax forest in the arid and semiarid soils of the Loess Plateau, China (Figure 1; Cheng, Han, & Kang, 2014), and the Q. wutaishansea species have cold tolerance and drought tolerance (Gao & Sun, 2005). Studies of vegetation variation and development have shown that, after herb and shrub communities, Pinus tabulaeformis, as a pioneer species, is generally the first established species and captures resources in disturbance openings (Fan, Guo, Wang, & Duan, 2014). The high recruitment rates of P. tabulaeformis are due to the high light intensity and temperatures of the bare land, and the seedlings are competitive, especially on the burnt ground (Gomez, 2004;Yang et al., 2018). Once P. tabulaeformis forests are formed, other, relatively shade-tolerant, species (i.e., Q. wutaishansea) are able to invade (Collins & Carson, 2004). The two species were found in a positive self-correlation and high intensity of competition (Li et al., 2008), which is consistent with those by Tilman (1994) andPedersen et al. (2012) who found that the competition intensity generally decreased with the growth of trees.
The opportunities for Q. wutaishansea to establish under the P. tabulaeformis were documented on the Loess Plateau, where a key strategy was establishment of shade-tolerant advanced reproduction (Rentch, Fajvan, & Rrjr, 2003). In the past 50 years, the P. tabulaeformis forest on the shady and semishady slope has been or is developing the Q. wutaishansea forest on the Loess Plateau. It has been suggested that pure P. tabulaeformis and Q. wutaishansea forests are at different stand development stages (Dang, Wang, & Wang, 2002), and the natural transformation of P. tabulaeformis F I G U R E 1 An organism photograph: this figure is for Quercus wutaishansea forests, which is the dominant or natural climax forest in the arid and semiarid soils of the Loess Plateau, China forests to mixed P. tabulaeformis and Q. wutaishansea forests has significant impacts on the underlying ecological interactions of the main species. Thus, a detailed investigation of forest spatial distributions may demonstrate and quantify the population mechanisms in natural secondary forests on the Loess Plateau. In this paper, we compare the composition, structure, spatial patterns, and interaction of dominant tree species among three different forest communities on the Loess Plateau. The objectives of this study were to determine how spatial patterns or aggregation change with scale and size class through the development process, and gain insight into tree coexistence and development trends. These results will provide the basis for future ecological restoration research and forest reconstruction, and offer a valuable reference for future planning of afforestation and ecological policies on the Loess Plateau and similar areas.

| Study area
Our study area was located on the southeastern Loess Plateau in northern Shaanxi Province. Detailed surveys were conducted on the Huanglong Mountains ( Figure 2). The entire forest zone shows extreme irregularity, a complicated geological structure, and fractured landforms, with the typical geographical structures of the Loess Plateau. The altitude ranges from 800 to 2,500 m, the frost-free period averages 175 days, annual precipitation is 611.8 mm (mostly from July to September), and the general climate is within the warm temperate zone (Guo, Li, He, Yang, & Sun, 2011). The steep and hilly topography mainly consists of granite and gneiss. The mean slope is 35°, and the mean soil depth is 50 cm. The soil is classified as mountain brown earth (Shi et al., 2004

| Data sets
Substituting space for time has been a widely used method for evaluating vegetation changes. It was initially used to distinguish successional stages in woody vegetation, and the method assumed that time was the only difference (Pickett, 1989 forests. Each forest was sampled on both a southern slope and a northern slope (with a due east-due west line as the boundary, the slope at the southern boundary is the sunny slope and the slope at the northern boundary is the shaded slope). With a canopy density of approximately 0.8 based on the adjacent gridding method, The location of the study area the plot was divided into 25 20 × 20 m small plots. In each small plot, the coordinates (x and y), diameter at breast height in 1.3 m (DBH), height (h), and the crown width of adult trees (DBH >3 cm) were measured (Table 1). In addition, 5 2 × 2 m 2 bush quadrats and 5 1 × 1 m 2 plant quadrats (1 × 1 m) were, respectively, arranged at the four corners and at the center of the subplot.

| Data analysis
Niche breadth is considered to be the determinants of species diversity and community structure, which was calculated using the Shannon index (Leibold, 1995;Magurran, 1988): where B i is the niche breadth of the ith species, with the range as [1/r, 1], r is the total number of resource states, and p ij is the proportion of the individuals of the ith species, which is associated with resource state j.
The importance value is a comprehensive quantitative indicator used to characterize the status and role of each species in the community, the larger the importance value of a species, the greater the dominance of the species in the plot, which was calculated with the following equation (Li et al., 2006): where RD is the relative dominance, RF is the relative dominance, RA is the relative dominance, n i is the number of individuals of the ith species, a i is the basal area at the height of 1.3 m belonging to the ith species, f i is the number of quadrats in which the ith species appeared, and S is the total number of species.
Ripley's K function K(r) is an important spatial pattern statistic, which is a cumulative distribution function within the distance of r (Getis & Franklin, 1987;Greig-Smith, 1952;Ripley, 1977). An alternative approach uses rings (or annulus) instead of circles, with a pair correlation function g(r) or the O-ring statistic O(r) (Wiegand, Moloney, Naves, & Knauer, 1999). The advantage of this alternative approach is that one can isolate specific distance classes, whereas the K(r) confounds effects at larger distances with effects at shorter distances (Condit et al., 2000;Stoyan & Penttinen, 2000).
In the present study, the univariate function g(r) was used to analyze the spatial distribution patterns in the categories of all trees and four size classes of trees. The function g(r) is related to Ripley's K function (Wiegand, Jeltsch, & Ward, 2000): where r is the distance (ring in g(r)) rad, A is the area of study plot, u ij is the distance between the focal tree (i) and its neighboring tree (j), n is the total number of points in the point pattern, u ij = 1, if u ij < r and 0 otherwise, and e ij is the weighting factor for eliminating edge effect correction.
For the univariate function g(r), the spatial scale was 0-50 m, and we used a ring width of one meter and used 199 Monte Carlo simulations of complete spatial randomness (CSR) to acquire pointwise critical envelopes, and the significance level was 0.01 (namely p ≤ 0.01). If the value of g(r) was above (or below) the upper limit of the confidence envelope, the spatial pattern indicates the clustered (or regular) pattern at a given distance r, and within the confidence intervals indicated random pattern.
To investigate the relationships between different life stages, we used bivariate pair correlation function g 12 (r), which is the extended g(r) function to multitype point patterns. g 12 (r) can be defined as the expected number of trees of life stage 2 at spatial scale r of an arbitrary tree of life stage 1, divided by the intensity of life stage 2 (Stoyan & Penttinen, 2000).
where u ij is the distance between the focal tree of pattern 1, and its neighboring tree of pattern 2, n 1 and n 2 are the total numbers of trees in the patterns 1 and 2, respectively (Wiegand & Molone, 2004).
For g 12 (r), we also used a ring width of one cell and used 199 Monte Carlo simulations of complete spatial randomness (CSR) to acquire pointwise critical envelopes and the significance level was 0.01 (namely p ≤ 0.01). The spatial associations between different growth stages were examined with the independent null model. If the value of g12(r) was above the upper (or below the lower) confidence limit, the relationship indicates that species are positively (or negatively) associated at the distance r, and within the confidence intervals indicated no interaction (Thioulouse, Chessel, Doledec, & Olivier, 1997;Wiegand et al., 2000).
We used package "spatstat" in R-3.4.2 to conduct all spatial analyses. In our analyses, a 1 m cell size was used, which preliminary analysis had suggested was a sufficient resolution to address our questions.

| Forest composition and structure
We identified 23 species, belonging to 12 families among the 911 live trees found in the 25 quadrats in P. tabulaeformis forests on the Loess Plateau. Twenty species, belonging to 14 families, respectively, among the 2,051 and 985 live trees were found in the 25 quadrats in Pinus + Quercus mixed forests and Q. wutaishansea forests (Tables 1   and 2). Twelve species were ubiquitous across all forests, and three of these were dominant species with importance values >5 ( The diameter class distribution showed that the numbers of individuals were gradually decreasing as the diameters increased ( Figure 3).
The diameter structure was relatively unitary and indicated that juvenile trees, ranging in diameter from 3.0 to 10.9 cm, accounted for trees in each stand. In P. tabulaeformis forest, the quantities of saplings for P. tabulaeformis and Q. wutaishansea were more than juveniles; and in Pinus + Quercus mixed forests and Q. wutaishansea forests, and sapling for Q. wutaishansea was more than juveniles (Figures 3 and 4). Dbh class distributions of P. tabulaeformis and Q. wutaishansea in each 1-ha plot indicated the natural regeneration ability of Q. wutaishansea was stronger than P. tabulaeformis, which was also related to the negative tolerance of Q. wutaishansea saplings.

| Spatial patterns at development stages
The spatial patterns among different development stages were estimated by the pair correlation function ( Figure 5). In the forest dominated by P. tabulaeformis, species distribution patterns transitioned from aggregated to random to uniform. Additionally, forests displayed an intensive aggregation pattern of small-scale at distances up to 14 m, a random pattern at distances from 15 to 19 m, and a uniform pattern above 19m. In Pinus + Quercus mixed forests, all trees were markedly aggregated at all distances (p < 0.01). In Q. wutaishansea forests, trees were aggregated at distances from 0 to 28 m and tended to display a random pattern at greater distances.

| Spatial patterns among life stages
We next analyzed the spatial distributions in different tree size classes ( Figure 6). In all three forests, juvenile trees showed an aggregated pattern at distances up to 10 m. In P. tabulaeformis forests and in Q. wutaishansea forests, juveniles displayed a random pattern at greater distances. In Pinus + Quercus mixed forests, however, juveniles showed an aggregated pattern up to distances of 13 m, and random between 29 and 36 m, and the value of g(r) was close to the upper broken fitting line.
The pattern of subadult trees in P. tabulaeformis forests and in Q. wutaishansea forests were aggregated at distances of 1-8 m; and at distances >8 m, distribution patterns were random. In Pinus + Quercus mixed forests, however, subadults were randomly distributed at all distances (p < 0.01).
Distribution of adult trees was random in all three forests over all distances. Saplings were randomly distributed in P. tabulaeformis forest and in Q. wutaishansea forests, but in Pinus + Quercus mixed forests displayed an aggregated pattern that declined with increasing distance.

| Spatial associations
Using bivariate statistics, we performed species associations among life stages in all forests at distances from 0 to 50 m. In P. tabulaeformis forests, our results suggested a positive association between juveniles and subadults from distances of 0-12 m, and no spatial associations between 23 and 50m (Figure 7).
Similarly, adults were positively spatially associated with juveniles and subadults at distances <1 m, but at other distances they were spatially independent. In Pinus + Quercus mixed forests, juveniles were negatively associated at low distances and random at high distances. Regardless of distance, spatial relationship showed all random patterns between adults and subadults. Juveniles and adults had a positive association across all distances. In Q. wutaishansea forests, juveniles showed random pattern with subadults at all distances. We observed few significant relationships between juveniles and adults. Subadults were spatially independent at distances <2 m and at the distances between 10 and 50 m. Between 3 and 6 m, subadults were negatively associated with adults.
F I G U R E 3 Dbh class distributions of all trees in each 1-ha plot. In this Figure, x-axes were displayed the four size classes, and the y-axes were displayed the percentage proportion of individual trees. (a) Pinus tabulaeformis forest, (b) Pinus + Quercus mixed forests, and (c) Quercus wutaishansea forest. The four size classes were "Juvenile" (3 < DBH < 10), "Subadult" (10 < DBH < 20), "Adult" (DBH ≥ 20), and "Sapling" is with heights <1.3 m was included

| Changes in population structure
It has been demonstrated that tree size was an important factor that affected stem growth and could be used to explain forest dynamics (Jang, Christopher, & Lim, 2013;Zhang, Gove, Liu, & Leak, 2001). In our research, the tree diameter class distribution curves across all three forests suggested that trees in all forests were highly dense as saplings, but that mortality was also high ( Figure 3). It has been shown that initial clustering structure for juveniles was an important factor affecting adult spatial structure (Zhao, Wang, Bai, Pan, & Wang, 2015), but environmental heterogeneity also was an important factor (Myster & Pickett, 1992;Schurr, Bossdorf, Milton, & Schumacher, 2004 Zou, Liu, and Wang (2002), Zou, Cheng, and Lei (1998) and Lei, Wang, Guo, and Zhu (2007).
It has been demonstrated that Q. wutaishansea is more shade tolerant than P. tabulaeformis (Zhu, 1999). This enabled Q. wutaishansea to persist in the shade of P. tablulaeformis, thereby forming more complex structures than when these species were growth separately. The population structure changed markedly between P. tabulaeformis forests and Q. wutaishansea forests, as expected.  (Fan et al., 2014). The net effect of these processes was a long-term reduction in P. tabulaeformis with consequences for the composition and F I G U R E 4 Dbh class distributions of Pinus tabulaeformis and Quercus wutaishansea in each 1-ha plot. In this Figure, x-axes were displayed the four size classes, and the y-axes were displayed the percentage proportion of individual trees. The four size classes were "Juvenile" (3 < DBH < 10), "Subadult" (10 < DBH < 20), "Adult" (DBH ≥ 20), and "Sapling" is with heights <1.3 m was included F I G U R E 5 Spatial patterns of all trees in the different successional stages using point pattern analysis method g(r): (a) Pinus tabulaeformis forest, (b) Pinus + Quercus mixed forests, (c) Quercus wutaishansea forest. We used 199 Monte Carlo simulations of complete spatial randomness (CSR) to obtain pointwise critical envelopes for g(r). We considered the pointwise Monte Carlo test significant at p ≤ 0.01. If g(r) within confidence intervals, then the spatial pattern at distance r is entirely random; if g(r) is above the upper confidence interval and below the lower confidence interval, then the spatial pattern indicates aggregated pattern and regular pattern, respectively (p < 0.01) structure of these forests (Gomez, 2004;Zou et al., 1998). There was a large-scale artificial P. tabulaeformis forest in the study area, and reasonable human intervention will gradually develop into the Pinus + Quercus mixed forests, and then develop into Q. wutaishansea forests (Wang, 1991).

| Spatial distribution
Our comparisons of spatial structure among different development stages suggested that stages differed greatly in spatial pattern.
Spatial patterns of pioneer species had less aggregation in older forests or in the less disturbed forests, as compared to forests that were younger or more disturbed. The action of competition had been shown to decrease aggregation and promoted shifting toward more regular or random patterns (López, Alcázar, & Ruiz, 2010;Martíneza et al., 2010).
Population distribution patterns depended significantly on spatial scale as has been reported elsewhere (Levin, 1992). We found that P. tabulaeformis forests exhibited small-scale aggregation and tended to form random patterns at large scales. This may be due to the lack of airbags, wings, and other flight aids; it was difficult to spread the seeds over long distances with the help of wind, and the seeds were concentrated around the mother plants. The complex terrain in the mountains also concentrates seeds in low places (Yang et al., 2018). This result was the same as previous studies, particularly other coniferous species which were aggregated in small scales and a randomly distribution in large scales (Cain & Shelton, 1995;Koukoulas & Blackburn, 2005;Shibata, Kikuchi, Tanaka, Sueyoshi, & Yoshimaru, 2009).
Under natural conditions, different tree sizes were typically found in separate canopy strata. Since stands may be subject to variable environmental conditions, this may affect distribution patterns (He, Legendre, & LaFrankie, 1997). We found that different stages of development had many juveniles which showed aggregated distributions at small scales. Aggregated spatial distributions are usually observed in naturally regenerating forests (Wang, Shao, & Shangguan, 2010), particularly in arid and semiarid forests (Xu, Xu, & Xie, 2012).
We found evidence of clumping in Pinus + Quercus mixed forests.
However, the degree of aggregation varied in forest types, and the distribution of aggregation could be explained by regeneration processes. Regeneration typically occurred near seed sources and was influenced by habitat heterogeneity (Jang et al., 2013).
At larger scales, the spatial patterns in these forests lost their strong aggregation and approached random for both juveniles and subadults. Unlike juveniles and subadults, random patterns were observed for adult trees at almost every scale, which may be because of intraspecific thinning due to high resource requirements and competition (Halpin & Lorimer, 2016;Tiphaine, Hugo, Frédérik, & Yves, 2014). Species that show aggregated patterns as young trees tend to present random patterns as the forest becomes more stable (Aldrich, Parker, Ward, & Michler, 2003). Our results indicated that density-dependent mortality of offspring was common. When dense plant communities are subject to intense competition, density-dependent mortality occurs in the process known as self-thinning (Yoda, Kira, Ogawa, & Hozumi, 1963). However, self-thinning, especially in natural forests, does not necessarily lead to random spatial patterns. For example, catastrophic natural events such as strong winds may result in random spatial patterns of adult P. tabulaeformis, which were more sensitive to wind damage due to shallow, wide root systems (Gao & Sun, 2005). However, such catastrophic events were rare and produce only local damage, so it was not easy to ensure that they would affect spatial patterns at a broad scale.
It was notable that the most limiting condition in the areas of semiarid and arid was water. Changes in forest cover (Whitford, 2003), soil surface cover, and texture affected water available to trees (Zou, Li, Xu, & Xu, 2010). For natural secondary forest, small-scale spatial patterns were mainly determined by biological characteristics such as species origin, initiation, and seed dissemination. In contrast, as spatial scale increased, the influence of topography, geomorphology, light, and moisture gradually replaced biological characteristics in importance, resulting in an alteration in spatial distribution patterns. These results lay the foundation for understanding the coexistence mechanism of major species in natural secondary forests and the impact of artificial replanting of forest stands. In the case of near-natural afforestation and management, it is necessary to pay attention to maintaining reasonable planting and growth density within and between species.

| Associations among life stage
Associations across life stages were mostly positive, and small-scale, probably attributed to prompt growth reactions to spatially limited light availability, which was the same opinions suggested by Hubbell, Ahumada, and Condit (2001). However, the positive association that we observed among adult trees and juvenile trees in P. tabulaeformis forests were somewhat surprising, as this pattern would be expected to vanish due to competition as trees grew to adult stage. Our results were in accordance with the findings of recent investigations (Martíneza et al., 2010). Additionally, Pacala (1997) suggested that intraspecific competition effects should be much more important than interspecific competition, which could promote species coexistence (Carson, & Pickett, 1990;Martíneza et al., 2010;Zhao, Kang, Guo, Yang, & Xu, 2012). It was interesting that the forests demonstrated self-thinning processes of intraspecific species, which could create favorable conditions for the invasion and settlement of other species. In order to promote the nature regeneration, the shrubs under the forest can be properly cleaned to provide more space for the sapling growth. At the same time, it is also possible to selectively harvest poorly growing trees and diseased trees to promote forest renewal. Moderate intervention on population distribution and adopting tending measures could promote forest succession.
Our analyses were based only on comparisons of spatial patterns and interspecific associations. However, it was noted that the 1-ha sample plots were representative of typical secondary natural forest composition, as compared to replicated forest series at smaller scales (Feroz, Yoshimura, & Hagihara, 2008;Whitfeld et al., 2014).
As there were similar factors in environmental heterogeneity among the three forests, the spatial pattern differences were likely caused F I G U R E 7 The bivariate statistic of the pair correlation function was used to analyze the spatial associations of the four size classes analyzed at the different growth stages at scales 0-50 m, which were Juvenile stage ( Our findings were somewhat expected, as interspecific differences in shade tolerance among tree species are key determinants of forest dynamics and structure, and the shade under the canopy of shade intolerant species is expected to facilitate establishment and growth of shade-tolerant species (Oliver & Larson, 1996). Therefore, young individuals of Q. wutaishansea would gradually dominate the natural space after a period of competition and elimination. Artificial disturbances, such as thinning or tending might accelerate this process, improve internal structure of forest and be conducive to forest development. Many factors influence the formation of forests, and no spatial pattern analysis can capture all of these factors. Therefore, to gain a better understanding of spatial pattern formation at different stages and to provide a framework for local forest management, the relationship among species (including P. tabulaeformis and B. platyphylla) and environmental factors should be further examined.

| CON CLUS IONS
Comparisons of spatial structure among our three forest types suggested that the most important influences on spatial patterns were the different forest development stages. Our study found that the differences in spatial patterns among various stages could inform mechanistic hypotheses explaining successional development. The seedlings of P. tabulaeformis were light -demanding and grew quickly with increasing light availability, while the seedlings of Q. wutaishansea were shade tolerant. Knowledge of this pattern can help to improve silvicultural practices. But due to the limited experimental conditions, the data in this study did not include trees with a DBH <3.0 cm. Therefore, discrepancies with other studies in the spatial distribution patterns and associations may be observed.
Stand development was difficult to evaluate because of its long timescale. Nevertheless, the length and structure of each stage in this study contributed to the site conditions, vegetation, and disturbance characteristics. Regenerated seedlings of Q. wutaishansea were found on slopes in P. tabulaeformis forests, and the proportion was more than that of P. tabulaeformis. In some stands, Q. wutaishansea occupied the second sublayer of the tree layer, demonstrating a notable tendency to replace P. tabulaeformis. In the process of forest management, artificial measures such as updating and arranging tree species can be adopted to accelerate the development of forests. It is worth noting that the mechanism of community stability structure needs to be further studied. The self-renewal of dominant species, the coexistence of species, and the community succession status of tree species needs to be long-term oriented observations.

ACK N OWLED G M ENTS
This project was supported by the National Natural Science

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R CO NTR I B UTI O N
Li Gu prepared figures and tables and wrote the main manuscript text; Kevin L. O'Hara was mainly responsible for the content of the thesis, revision, and language adaptation; Weizhong Li was responsible for field survey and data collection; Zhiwen Gong was the corresponding author. And all authors reviewed the manuscript.

DATA ACCE SS I B I LIT Y
Data available from the Dryad Digital Repository: https ://doi. org/10.5061/dryad.hk31vm7. Please see lines 397.