Using Count Data Models to Predict Epiphytic Bryophyte Recruitment in Schima superba Gardn. et Champ. Plantations in Urban Forests

Epiphytic bryophytes are known to perform essential ecosystem functions, but their sensitivity to environmental quality and change makes their survival and development vulnerable to global changes, especially habitat loss in urban environments. Fortunately, extensive urban tree planting programs worldwide have had a positive effect on the colonization and development of epiphytic bryophytes. However, how epiphytic bryophytes occur and grow on planted trees remain poorly known, especially in urban environments. In the present study, we surveyed the distribution of epiphytic bryophytes on tree trunks in a Schima superba Gardn. et Champ. urban plantation and then developed count data models, including tree characteristics, stand characteristics, human disturbance, terrain factors, and microclimate to predict the drivers on epiphytic bryophyte recruitment. Different counting models (Poisson, Negative binomial, Zero-inflated Poisson, Zero-inflated negative binomial, Hurdle-Poisson, Hurdle-negative binomial) were compared for a data analysis to account for the zero-inflated data structure. Our results show that (i) the shaded side and base of tree trunks were the preferred locations for bryophytes to colonize in urban plantations, (ii) both hurdle models performed well in modeling epiphytic bryophyte recruitment, and (iii) both hurdle models showed that the tree height, diameter at breast height (DBH), leaf area index (LAI), and altitude (ALT) promoted the occurrence of epiphytic bryophytes, but the height under branch and interference intensity of human activities opposed the occurrence of epiphytic bryophytes. Specifically, DBH and LAI had positive effects on the species richness recruitment count; similarly, DBH and ALT had positive effects on the abundance recruitment count, but slope had a negative effect. To promote the occurrence and growth of epiphytic bryophytes in urban tree planting programs, we suggest that managers regulate suitable habitats by cultivating and protecting large trees, promoting canopy closure, and controlling human disturbance.


Introduction
Epiphytic or corticolous bryophytes grow on the bark of living trees and shrubs [1] and are found to be widely distributed from the northern forests in the northern hemisphere to the temperate forests in the southern hemisphere, including tropical forests. Recent research has shown that epiphytic bryophytes play vital roles in forest ecosystems, such as storing atmospheric water and maintaining the microhydrological cycle in forests [2]; nutrient cycling in terrestrial ecosystems through atmospheric sedimentation, rainfall, and symbiosis with nitrogen fixing organisms [3]; providing important habitats for protozoa, invertebrates, arthropods, and other small creatures [4]; and as bio-indicators for environmental monitoring because of their special physiological structure [5,6].
The colonization and subsequent development of epiphytic bryophyte diversity are regulated by factors at the tree, stand, and global scales [7,8]. At the tree level, many studies have emphasized that the total tree height (H) [9], diameter at breast height (DBH) [10], canopy openness [11], age [12], and bark characteristics [13] are the driving factors of the diversity of epiphytic bryophytes. At the stand level, diversified microhabitats drive the development of epiphyte diversity that results from stand characteristics such as tree species composition [14], the continuity of the forest area [15], forest age [16], and tree density [17]. Epiphytic bryophytes are also affected by global changes caused by human activities, such as climate change [18], rising atmospheric carbon dioxide levels [19], nitrogen deposition [20], sulfur dioxide pollution [21]. In addition, most studies are conducted in native forest ecosystems, and relatively few studies involve urban forest ecosystems.
Urban environments have been shown to have important effects on epiphytic bryophytes, e.g., air pollution [22], urban heat island effect [23], habitat fragmentation [24], and other problems associated with urbanization can affect bryophyte growth and development. Most previous studies have focused on surveying epiphyte diversity [25], determining floristic changes of epiphytic bryophytes [26], and evaluating atmospheric conditions in urban areas [5,27]. Presently, an increasing number of countries are committed to promoting urban tree planting [28][29][30], which could provide potential habitats for epiphytic bryophytes, especially trees in parks. However, there is lack of knowledge about the recruitment of epiphytic bryophytes on urban trees.
In previous studies, numerous statistical techniques were used to assess the growth and development of epiphytic bryophytes, such as redundancy analysis and canonical correspondence analysis [7], the generalized linear model [31], structural equation models [32], linear mixed effect models [33], multiple regression models [34], and diverse combinations of multiple methods [14]. However, few studies have paid attention to the non-occurrence of epiphytic bryophytes. Therefore, we used count data models to solve this problem. Count data models are used to study the relationship between the explanatory variables and response variables. In count data models, the number of occurrences of bryophytes is considered as the response variable. Numerous studies have used a two-stage approach to solve the problem of a large number of zero values [35][36][37]: the first step is to calculate the probability of occurrence through a logistic model, and the second step is to predict differences in occurrence of the count by Poisson or negative binomial model. Presently, count data models are applied in medical science [38,39], biological statistics [35,40], disaster predictions [41,42] and other fields, but few studies for epiphytic bryophytes in other ecosystem [14]. To our knowledge, this study is the first to discuss the recruitment of epiphytic bryophytes in an urban forest plantation in a subtropical area using count data models.
In this study, we aimed to understand the occurrence and growth of epiphytic bryophytes and the possible driving factors after afforestation in subtropical urban forests. We assessed the epiphytic bryophytes on S. superba trees planted at the beginning of the 21st century in the Yangtaishan Forest Park in Shenzhen, China. The objectives of our study were: (i) to investigate the colonization of bryophytes on tree trunks; (ii) to document the factors that influence species richness and abundance of epiphytic bryophytes on trees; and (iii) to develop and compare Poisson, Negative binomial (NB), Zero-inflated Poisson (ZIP), Zero-inflated negative binomial (ZINB), Hurdle-Poisson (HP), and Hurdle-negative binomial (HNB) models to predict epiphyte recruitment.

Study Area
The study area is located in the Yangtaishan Forest Park (22 • 38 -40 N and 113 • 55 -59 W, altitude: 33.4-587.1 m a.s.l.), which covers a total area of 2850 ha in the west of Shenzhen City, in Guangdong province, southern China ( Figure A1). The climate in this area is a subtropical marine climate characterized by hot and rainy summers. According to the climate data from 1981 to 2010 [43], the average monthly temperature was highest in July (28.9°C) and lowest in January (15.4°C), the average monthly rainfall was highest in August (354.4 mm) and was smallest in January (26.4 mm), and the average monthly relative humidity was highest in June (80%) while lowest in December (64%). Throughout the study area, the dominant soil type is latosol and the dominant vegetation type is the evergreen broad-leaf mixed forest [44].
To eliminate the effect of tree species diversity on epiphytic bryophytes recruitment, we investigated an S. superba forest belt with an altitude gradient from 153.5 m to 534.0 m through the east and west main entrances of the park. The S. superba forests were planted from 1999 to 2005 for purposes of forest fire prevention. The trees were initially planted at a spacing of 2 × 2 m.

Site Selection
We selected 21 locations 300 to 500 m apart [32], along the S. superba forest belt ( Figure 1a). Two 10 × 20 m sample plots were established on both sides of the hiking trail at each location. In each plot, we randomly selected 30 healthy trees excluding any dead or leaning trees and only chose trees with DBH greater than 10 cm [45]. A total of 1260 S. superba trees were investigated in 42 plots. On each tree, we established 16 10 × 10 cm subplots in all: one subplot per orientation (East, South, West, and North) at each height (0.3, 1.1, 1.5, and 1.8 m up from the base of the tree) ( Figure 1c). In addition, we used the park square as a microclimate control location to represent urban environment conditions outside forest.

Study Area
The study area is located in the Yangtaishan Forest Park (22°38′-40′ N and 113°55′-59′ W, altitude: 33.4-587.1 m a.s.l.), which covers a total area of 2850 ha in the west of Shenzhen City, in Guangdong province, southern China ( Figure A1). The climate in this area is a subtropical marine climate characterized by hot and rainy summers. According to the climate data from 1981 to 2010 [43], the average monthly temperature was highest in July (28.9 ℃) and lowest in January (15.4 ℃), the average monthly rainfall was highest in August (354.4 mm) and was smallest in January (26.4 mm), and the average monthly relative humidity was highest in June (80%) while lowest in December (64%). Throughout the study area, the dominant soil type is latosol and the dominant vegetation type is the evergreen broad-leaf mixed forest [44].
To eliminate the effect of tree species diversity on epiphytic bryophytes recruitment, we investigated an S. superba forest belt with an altitude gradient from 153.5 m to 534.0 m through the east and west main entrances of the park. The S. superba forests were planted from 1999 to 2005 for purposes of forest fire prevention. The trees were initially planted at a spacing of 2 × 2 m.

Site Selection
We selected 21 locations 300 to 500 m apart [32], along the S. superba forest belt ( Figure 1a). Two 10 × 20 m sample plots were established on both sides of the hiking trail at each location. In each plot, we randomly selected 30 healthy trees excluding any dead or leaning trees and only chose trees with DBH greater than 10 cm [45]. A total of 1260 S. superba trees were investigated in 42 plots. On each tree, we established 16 10 × 10 cm subplots in all: one subplot per orientation (East, South, West, and North) at each height (0.3, 1.1, 1.5, and 1.8 m up from the base of the tree) ( Figure 1c). In addition, we used the park square as a microclimate control location to represent urban environment conditions outside forest.

Bryophytes Data Collection
According to the technical standards of the Ministry of Ecology and Environment [46], we investigated the number of bryophyte species and the coverage of each bryophyte in each subplot. Then, we determined the total number of bryophyte species in the 16 subplots on each selected tree

Bryophytes Data Collection
According to the technical standards of the Ministry of Ecology and Environment [46], we investigated the number of bryophyte species and the coverage of each bryophyte in each subplot. Then, we determined the total number of bryophyte species in the 16 subplots on each selected tree to express species richness (the response variable in the models). Specimens of bryophyte species were transported to the laboratory for storage and further identification. The bryophyte reference works of Wu and Zhang [47] and Zhang [48] were used as the main references for species identification. Coverage of each bryophyte per subplot was registered as the number of small squares by using a square grid with 100 squares (Figure 1d). Total coverage of each species was calculated in the 16 subplots on each tree, and then divided by the total coverage of the 16 subplots (1600 cm 2 ) as the abundance. To better use the count data models, the abundance on each tree was classified into grades: 0 (0%), 1 (0%-1%), 2 (1%-5%), 3 (5%-10%), 4 (10%-20%), 5 (20%-30%), 6 (30%-40%), etc. [49,50]. The abundance grade was chosen as the response variable in the models. The survey of bryophytes was performed between November 2018 and March 2019.

Environmental Data Collection
Twelve environmental factors were investigated to reflect tree characteristics, stand characteristics, terrain factors, microclimate, and human effects (Table 1). H, DBH, crown width (CW), and height under the lowest branch (HB) were measured for all trees. HB, defined as the vertical height from the ground to the lowest point of the crown branch, was measured using a TruPulse 200 Laser Rangefinder (Laser Technology, Centennial, CO, USA), as was H. DBH, defined as the trunk diameter at 1.3 m above the ground surface, was measured with a ruler. CW was measured using the vertical projection of tree crown on the ground. Leaf area index (LAI) and crown density (CD) were measured using an LAI-2200C plant canopy analyzer (Li-Cor, Lincoln, NE, USA). Terrain factors included altitude (ALT), slope (SLP), and northness index (NI). The aspect was converted to NI by using the cosine transformation [51]. The microclimate of each sample was recorded at 1 h intervals for 6 days without rain in November 2018, March 2019, and July 2019. The air temperature and air relative humidity were measured 1.8 m above the ground using a TR-74Ui-S temperature and humidity data logger (T&D Corporation, Matsumoto, Japan). Since we could not measure one day climate data simultaneously for all plots, we used relative climate data to reflect the climatic condition. Air temperature and humidity in each plot, divided by that in the control locality, were regarded as the relative air temperature (RT) and the relative air relative humidity (RH), respectively. Therefore, RT and RH had no unit. In addition, we counted the number of people entering the woodlands in each plot for 6 days in November 2018, March 2019, and July 2019. Then, we divided this number by the number of people with the maximum count plot and considered it as a proxy of interference intensity of human activities (IHA).

Data Analysis
To explore the tree-trunk colonization preference of bryophytes, we investigated the distribution of bryophytes at the subplot level. We calculated the total number of bryophyte species and coverage per subplot. Then, we divided this by the total number of 16 subplots and considered this as the percentage of species richness and abundance per subplot. We analyzed the preference of bryophytes for different height and direction of the trunk by comparing the values of proportion. The calculations and graphs were completed in Microsoft Excel 2016 (Microsoft, Redmond, WA, USA).
To model epiphytic bryophyte recruitment, we developed and compared different common count data models because all response variables (species richness and abundance) were in count form. The ZIP, ZINB, HP, and HNB models used a two-stage approach to accurately model recruitment instead of the ordinary least-squares methods [52]. First, a logistic function was used to predict whether epiphytic bryophytes occur on trees; from the model results, the factors that affect the probability of epiphyte colonization on trees can be determined. Then, the Poisson or NB models were used to describe the positive counts using the least squares method [35]. On the basis of the occurrence of epiphytic bryophytes, this step determined the growth of the recruitment count. Before modeling, ALT and LAI were log-transformed.
Modeling involved the following phases: (1) Selection of variables: we analyzed the multicollinearity between explanatory variables using a variance-inflation factor (VIF) test to determine the most representative explanatory variables and avoid potential issues with multicollinearity in the modeling process. Based on the test results, variables with VIF values higher than 5 were considered for deletion from the model [35,53]. (2) Model fitting: at the tree level, six model types were generated to model the species richness and abundance of epiphytic bryophytes with the remaining variables. To simplify the models, we progressively removed the least important variables until only the significant variables remained [54]. (3) Model selection: in our models, candidate explanatory variables were stepwise introduced into the models [55]. We compared model performance using the Akaike information criterion (AIC). In this process, the lowest AIC value indicated the most suitable model [56]. Vuong tests were performed for the non-nested fitted models with similar AIC values [57]. (4) Model goodness-of-fit: we used R-squared and residual diagnostics to evaluate the goodness-of-fit of each model. We calculated the coefficient of determination between predicted and measured values and chose to use the Pearson residuals because the original residuals were heteroscedastic and asymmetric for count data [58]. Modeling was established using packages car and pscl of the statistical software R version 3.6.1 (R Foundation for Statistical Computing, Vienna, Austria) [59].

Overall Species Diversity and Distribution
Overall, 17 different bryophyte species (10 mosses, 7 liverworts) were found on the trunks of S. superba trees. From the overall survey data, bryophytes occurred on 599 of 1260 sampled trees. The frequency of occurrence of Sematophyllum subhumile (Müll. Hal.) M. Fleisch. was the highest of all bryophytes registered (Table A1). The total number of epiphytic bryophyte species per tree was: 0 (661 trees), 1 (355 trees), 2 (196 trees), 3 (39 trees), and 4 (9 trees), i.e., a five-point scale classification system (Figure 2a). The abundance grade of epiphytic bryophytes was: 0 (661 trees), 1 (286 trees), 2 (210 trees), 3 (57 trees), 4 (31 trees), 5 (8 trees), 6 (7 trees), i.e., a seven-point scale classification system (Figure 2b). The bryophyte species, their tree and subplot level frequency values, and their total coverage are listed in Table A1.  Table A1. As can be seen in Figure 3, there were differences in bryophyte colonization in terms of the direction and height on the trunk. Regarding the direction, the percentage of the bryophyte species richness and coverage both reached maximum values on the side of the trunk facing north (88.24% and 39.19%, respectively). Regarding the trunk height, the percentage of the species richness and coverage both reached maximum values at 0.3 m above the ground (100% and 59.22%, respectively) ( Figure 3). Among the 16 subplots on the trunk, maximum species richness and coverage both occurred the subplot at the base of the trunks (0.3 m above the ground) in the north-facing direction.  As can be seen in Figure 3, there were differences in bryophyte colonization in terms of the direction and height on the trunk. Regarding the direction, the percentage of the bryophyte species richness and coverage both reached maximum values on the side of the trunk facing north (88.24% and 39.19%, respectively). Regarding the trunk height, the percentage of the species richness and coverage both reached maximum values at 0.3 m above the ground (100% and 59.22%, respectively) ( Figure 3). Among the 16 subplots on the trunk, maximum species richness and coverage both occurred the subplot at the base of the trunks (0.3 m above the ground) in the north-facing direction.  Table A1. As can be seen in Figure 3, there were differences in bryophyte colonization in terms of the direction and height on the trunk. Regarding the direction, the percentage of the bryophyte species richness and coverage both reached maximum values on the side of the trunk facing north (88.24% and 39.19%, respectively). Regarding the trunk height, the percentage of the species richness and coverage both reached maximum values at 0.3 m above the ground (100% and 59.22%, respectively) ( Figure 3). Among the 16 subplots on the trunk, maximum species richness and coverage both occurred the subplot at the base of the trunks (0.3 m above the ground) in the north-facing direction.

Species Richness Recruitment Models
Based on the results of VIF (Table A2), CD, RT, and RH were deleted from species richness and abundance recruitment models. As can be seen in Table 2, the HP and HNB models were more suitable than the other models because the AIC values for these models were both smaller than those for the Poisson, NB, ZIP, or ZINB models [60]. The Vuong test also indicated that the HP and HNB models were both more suitable than the ZIP and ZINB models. Furthermore, we have no strong evidence to suggest any significant difference between the HP and HNB models, with the results of both models being highly consistent. The goodness-of-fit of the two models was tested by Pearson residuals. The residuals both ranged from −2 to 3 for the HP and HNB models and the distribution of the residuals was comparatively uniform (Figure 4a,b). Finally, the R-squared values of the HP and HNB models were both 0.43. Values in parentheses represent standard deviation. The asterisks indicate significant difference: * p < 0.05, ** p < 0.01, *** p < 0.001.
In the zero component of the hurdle models, H, DBH, LAI, and ALT were found to promote the occurrence of bryophytes on S. superba trees, whereas HB and IHA had a negative influence. All factors were significantly correlated at the 0.001 level of significance. In the positive count component of the hurdle models, only DBH and LAI were found to have a positive effect on species richness recruitment count at the 0.001 and 0.05 levels of significance, respectively.

Abundance Recruitment Models
The results of modeling abundance of epiphytic bryophytes and explanatory variables are shown in Table 3. As seen in the models, the AIC values and Vuong test for species richness indicated that the hurdle models (HP and HNB) were more suitable than any other model. The Pearson residual plots, which showed similar results to the species richness residuals, showed that, except for one outlier, most of the residuals of the HP and HNB models ranged from −2 to 3 and the distribution of the residuals was comparatively uniform (Figure 4c and d). The R-squared values of the predicted values and the actual values of the two models were both 0.47.  In the zero component of the hurdle models, H, DBH, LAI, and ALT were found to promote the occurrence of bryophytes on S. superba trees, whereas HB and IHA had a negative influence. All factors were significantly correlated at the 0.001 level of significance. In the positive count component of the hurdle models, only DBH and LAI were found to have a positive effect on species richness recruitment count at the 0.001 and 0.05 levels of significance, respectively.

Abundance Recruitment Models
The results of modeling abundance of epiphytic bryophytes and explanatory variables are shown in Table 3. As seen in the models, the AIC values and Vuong test for species richness indicated that the hurdle models (HP and HNB) were more suitable than any other model. The Pearson residual plots, which showed similar results to the species richness residuals, showed that, except for one outlier, most of the residuals of the HP and HNB models ranged from −2 to 3 and the distribution of the residuals was comparatively uniform (Figure 4c,d). The R-squared values of the predicted values and the actual values of the two models were both 0.47.
The results of the hurdle models were the same as those from the zero component in the species richness recruitment model because the same data were used in the zero component of both the species richness and abundance recruitment models. In the positive count component, only DBH and LAI were found to have a positive effect on the abundance recruitment count, while the opposite relationship was observed between SLP and the abundance recruitment count. Table 3. Parameter estimates and fit statistics of Poisson, Negative binomial (NB), Zero-inflated Poisson (ZIP), Zero-inflated negative binomial (ZINB), Hurdle-Poisson (HP), and Hurdle negative binomial (HNB) models for abundance of epiphytic bryophytes.

Overall Species Diversity and Distribution
The number of epiphytic bryophyte species in the S. superba plantations was almost half of the number found in the Yangtaishan Forest Park (38 species, unpublished data). In addition, eight bryophytes were found on 10 S. superba trees with an average coverage of 3.98% in Heishiding Nature Reserve [61]. In the present study, we found 17 bryophytes on 1260 S. superba trees with an average coverage of 1.51%. These findings demonstrate that plantations in urban environments in subtropical areas can provide habitats for bryophytes.
In the present study, assessment of the distribution of bryophytes on parts of the trunks facing different directions revealed that epiphytic bryophytes mostly preferred to colonize the north side (i.e., shaded side) of the trunks. Although some studies have indicated that specific species might have a preference for one side of the trunk [62], in general, most species prefer the shaded side [63]. The average photosynthetic saturation of bryophytes is lower than that of vascular plants; therefore, the lower light intensity on the shaded side may provide a microhabitat that is more conducive to growth [64]. In addition, such microhabitats may enable bryophytes to avoid desiccation because of the reduced water evaporation [65].
Consistently with the aforementioned reports, the base of the trunks (0.3 m) had both the highest species number and coverage of epiphytic bryophytes. This colonization preference may be attributed to the relatively higher RH, lower RT, and more stable conditions in the microclimate at the base of trunk [66]. Thus, the bryophyte colonization preference observed in this study was consistent with findings reported by similar studies in temperate and subtropical forests [6,67].
Considering the distribution of bryophytes on parts of the trunks facing different directions and at different heights, our results indicate that the position of the trunk with moist and stable microclimates may be the preferred place for bryophytes in urban plantations in subtropical areas.

Factors Affecting Epiphytic Bryophyte Recruitment
One of the purposes of this study was to predict the probability of epiphytic bryophyte recruitment on trees; the zero part of the hurdle models explained this issue. H, DBH, HB, LAI, ALT, and IHA were found to be the most important factors determining the colonization success rate of epiphytic bryophytes. H determined the microclimate changes along the whole vertical gradient, to attract the colonization by epiphytic bryophytes with different microhabitat preferences [9]. DBH had a positive effect on the occurrence of epiphytic bryophyte, consistently with the findings reported by numerous previous studies [9,17,68]. This may be due to the fact that trees with larger DBH values have greater bark moisture content and enriched bark texture [69]. Interestingly, the probability of occurrence of epiphytic bryophytes decreased with increasing HB; likely because the higher the HB, the lighter the exposure and the lower the RH on the trunk [67]. LAI and ALT both affected the probability of epiphytic bryophyte colonization by regulating the microclimate of forests, which was confirmed by the correlations among LAI, ALT, RT, and RH (Table A3). As each factor directly or indirectly affects relative humidity, epiphytic bryophytes may have a preference for cool and humid microclimate. Although bryophytes have poikolohydric characteristics, which enable them to adapt to drought conditions, atmospheric moisture content was still found to be the most important factor that limiting their growth because of their lack of roots and effective water conducting tissues. Therefore, bryophytes, especially in the early stages of colonization, preferred the cool and moist microclimate prevalent on the lower, north-facing parts of S. superba trunks [7,31]. In addition, we found that human disturbance generated an inhibitory effect on the colonization of epiphytic bryophytes, which consistent with other research [70], human disturbances seeming affect microclimate stability.
A positive effect of DBH on recruitment count was found in the positive count components of both the models of species richness and abundance recruitment. DBH become the only common factor for the occurrence and the growth of species richness and abundance of epiphytic bryophytes. As discussed previously, an increase in DBH leads to an increase in bark surface area and a richer bark texture, thereby helping to provide more diversified habitats to attract colonization by more epiphytic bryophytes species [71]. Large DBH implies that trees are old enough to have provided a longer time for colonization by more epiphytic bryophyte species [72]. In addition, the moisture content of the bark may increase with increasing DBH [10], which is also conducive to the growth of epiphytic bryophytes. Another positive factor on recruitment count for species richness was LAI. LAI is strongly associated to CD (Table A3), which might lessen light exposure and increase RH in the forest, and thus promoting colonization by more species. These findings were consistent with a study conducted in a dry forest, whereby the species density of epiphytic bryophytes was highest in areas with a closed canopy [73].
In addition, ALT and SLP had different effects on bryophyte abundance. ALT had a positive effect on abundance recruitment count, but not on species richness. That was probably caused by the mid-altitude bulge phenomenon, which species richness of bryophytes peaked at a mid-altitude [74,75]. Conversely, SLP had a negative effect on bryophyte abundance, presumably because light exposure within the forest increases with increasing SLP, which adversely affects the growth of bryophytes [76]. On the other hand, we found that the correlation between SLP and IHA was 0.34 (p < 0.01) (Table A3). It is possible that disturbances caused by people walking in the nearby S. superba forest instead of the stone pathways might have led to the instability of the microclimate, thus affecting the growth of bryophytes. Previous studies have confirmed the importance of aspect for the growth of bryophytes [77]; however, the experimental plots in this study were located near the ridge line. In such a case, as the difference in light exposure and microclimate between north and south orientations was likely negligible, aspect had little effect on the growth of bryophytes.

Comparison Among Basic Recruitment Models
According to the results of the AIC and Vuong tests, our study suggests that the hurdle models perform better than other models in predicting epiphytic bryophyte recruitment. Compared with the Poisson and NB models, the hurdle models yielded better simulation results because the hurdle models divided the bryophyte recruitment process in two phases: occurrence and growth. Compared with zero-inflated models, the hurdle models performed better because the zero component of the hurdle models came only from the data without the occurrence of epiphytic bryophytes (the structure zero part) [78]. However, the zero component of the zero-inflated models may have come from both the structure zero portion, and from the positive count component in the Poisson or negative binomial distribution [79]. Thus, by dealing only with the zero counts for epiphytic bryophyte recruitment, our study found that the hurdle models were slightly more suitable than other models.
In this study, the AIC and Vuong tests did not provide sufficient evidence to determine any difference between the HP and HNB models. The Poisson model requires the variance of the outcome, which is assumed to equal its mean, while the NB model is better for overdispersion data [80]. Thus, for example, some tree recruitment studies showed that the negative binomial model was better, probably because the data was exceedingly discrete [35,53], whereas, in the present study, species richness and abundance data showed no overdispersion and the mean and the variance for species richness (0.72 and 0.89, respectively) and abundance (0.86 and 1.15, respectively) were relatively close (Table 1 and Figure 2). Consistently, the goodness-of-fit of the HP and HNB models, verified by Pearson residuals and R-squared values, did not show any difference. Therefore, our results show that the hurdle models (HP and HNB) both had good fitting abilities for continuous counting data with little difference between the mean and the variance.

Conclusions and Suggestions
To our knowledge, this is the first study to clarify the colonization patterns and influential factors of epiphytic bryophytes using count data models. Our results indicate that epiphytic bryophytes in urban environments prefer to colonize the shaded, humid parts of the trunk. Among the six count data models, both hurdle models, HP and HNB, predicted epiphytic bryophyte recruitment more suitably. These models showed that H, DBH, LAI, ALT, HB, and IHA were the main variables determining the occurrence of epiphytic bryophytes on trees of S. superba. Furthermore, DBH and LAI were found to positively affect the species richness recruitment count of epiphytic bryophytes, while DBH, ALT, and SLP were shown to play important roles in influencing the abundance recruitment count.
Our findings indicate that greater epiphytic bryophyte recruitment is more likely in habitats with large trees; cool, moist, and stable microclimates; and less human disturbance. Therefore, these factors are particularly important for urban managers to consider. To protect and restore the diversity of epiphytic bryophytes in urban forests, cultivating and protecting large trees, promoting canopy closure, controlling human disturbance, and other measures aimed to enhance suitable habitats for epiphytic bryophytes are needed in the early stages of afforestation efforts. If we wish to restore biodiversity and promote successful ecological succession of epiphytic bryophytes in urban areas, the factors that promote epiphyte recruitment should be considered. In addition, other potential factors should be thoroughly investigated, and continuous monitoring should be conducted in urban forests to further our understanding of these expanding environments.       Asterisks indicate significant difference: * p < 0.05, ** p < 0.01.