The mid‐domain effect and habitat complexity applied to elevational gradients: Moss species richness in a temperate semihumid monsoon climate mountain of China

Abstract The utility of elevational gradients as tools to test either ecological hypotheses and delineate elevation‐associated environmental factors that explain the species diversity patterns is critical for moss species conservation. We examined the elevational patterns of species richness and evaluated the effects of spatial and environmental factors on moss species predicted a priori by alternative hypotheses, including mid‐domain effect (MDE), habitat complexity, energy, and environment proposed to explain the variation of diversity. Last, we assessed the contribution of elevation toward explaining the heterogeneity among sampling sites. We observed the hump‐shaped distribution pattern of species richness along elevational gradient. The MDE and the habitat complexity hypothesis were supported with MDE being the primary driver for richness patterns, whereas little support was found for the energy and the environmental factors.

Numerous hypotheses have been proposed to explain species richness patterns along elevational gradients, but no one is consistently supported with empirical data (Nascimbene & Marini, 2015;Raabe et al., 2010;Spitale, 2016). Specifically, the mid-domain effect (MDE) indicates that if species' ranges are distributed randomly within a bounded domain, more ranges will overlap in the middle of the domain than at the edges which will produce a humpshaped pattern of species richness (Colwell & Hurtt, 1994;Colwell & Lees, 2000). The energy hypothesis proposes that higher ambient energy and productivity often results in higher species diversity (Hawkins et al., 2003). The environment hypothesis proposes that species richness patterns are generated by the climatic factors such as rainfall, temperature, and water availability (Heaney, 2001;McCain, 2007;Sánchez-Cordero, 2001). Moreover, habitat complexity has also been regarded as a potential driver of species richness (Brown, 2001;Wu et al., 2013).
Although the mechanism of geographic variation in species richness is important and has been explored by ecologists for decades, there are still limitations of elevational richness patterns. Species are generally not homogeneously distributed along elevational gradients, and the heterogeneity in biodiversity within (α) and among (β) sampling sites cannot be revealed by the elevational richness pattern. As envisioned by the combination of additive diversity partitioning and species-area relationship, β-diversity among sampling sites may partly be explained by a factor gradient (Gao & Perry, 2016;Golodets et al., 2011;Zajac et al., 2013); thereby, we suggest elevational richness pattern alone cannot quantify how much of the total β-diversity is due to elevation (β elevation ) and how much is due to other factors (β replace ). Moreover, the comparison of the diversity within (α) and among (β) sampling sites and the contributions made by elevation (β elevation ) and other factors (β replace ) are important for strategic conservation planning. A low α-diversity with a high β-diversity suggests that species assemblages are heterogeneous and species are often specific to individual sampling sites, while a high α-diversity with a low β-diversity indicates that species assemblages are homogenous and species within each sampling site are a subsample of the same species pool (Francisco-Ramos & Arias-González, 2013). A high β elevation with a low β replace indicates that species richness varies in a more predictable manner determined by factors that have a strong association with elevation, while a low β elevation with a high β replace suggests that factors such as speciation, dispersal, and extinction have a greater role in influencing patterns of β-diversity (Rahbek, Borregaard, Antonelli, et al., 2019).
In contrast to many vascular plants, mosses are dispersed by means of small spores and establish new populations in distant localities.
They colonize almost all terrestrial habitats, exhibit less frequent speciation, and have a long evolutionary history. Because of these unique features, our aims were to (a) depict the species richness pattern of mosses along elevational gradient, (b) evaluate the importance of four ecological hypotheses in predicting variation of moss species diversity along elevational gradient, and (c) examine how much contribution that elevation made toward explaining the among-sampling heterogeneity.

| Study area
This study was conducted on the Mt Tuofeng (maximum elevation 2,282 m a.s.l.). It sits within the Tuoliang National Reserve, in the middle of the Taihang Mountains in central west Hebei Province, China (38°33′-38°45′N, 113°41′-113°53′E). The study area is in the transition between warm and cold temperate zone and generally has a semihumid semiarid continental monsoon mountain climate with four distinct seasons, abundant sunshine, large temperature difference between day and night, moderate rainfall, and an annual average temperature (AT) of 8.0°C.

| Sampling and species identification
The present study was conducted along the elevational gradient of the Mt Tuofeng between 923 and 2,282 m a.s.l. during July-September 2018 within Tuoliang National Reserve. We randomly selected 73 sampling sites (10 × 10 m) from three transects along the elevational gradient to cover all types of vegetation with an equal elevational distance (c. 19 m) and over 100 m apart between each other. In each site, we collected all moss species from the ground to two meters above the ground. Moss specimens were taken back to the laboratory of Hebei Normal University where all species were identified from October 2018 to May 2019. Finally, the outcome of species occurrence for each sampling site is summarized in Table S1.

| Ecological variables: MDE, habitat, and climate
We utilized the MDE model in RangeModel ver. 5 (Colwell, 2006) to test the MDE. We employed the discrete domain analysis for the sampling sites, which in our study were discrete and evenly spaced.
Species richness data for each sampling site were compared with null model predictions using a Monte Carlo simulation of species richness curves to evaluate the explanatory power of the MDE on the species richness pattern. Simulated curves were based on empirical range sizes within a bounded domain, using the analytical stochastic models of Colwell and Hurtt (1994;Colwell & Lees, 2000).
We conducted 100,000 Monte Carlo simulations of empirical range sizes sampled without replacement (i.e., the randomization procedure) to calculate the mean expected species richness and their 95% confidence intervals for each sampling site.
We applied two habitat indices, including vegetation type (VT) and community type (CT; Chen, 1958Chen, , 1963, and counted the total number of VT and CT to quantify habitat complexity for each sampling site. We in total categorized six VTs, including coniferous forest, broad-leaved forest, bush wood, shrub grass, grass, and meadow, and four CTs, including aquatic community, stone community, soil community, and woody community. We recorded the VT for each sampling site and the CT occupied by each moss specimen. We obtained climate data from the WorldClim v2 database (Fick & Hijmans, 2017) in 30-arc-second (c. 1 km 2 ) digital maps, including solar radiation (SR), annual precipitation (AP), annual AT, and wind speed (WS). The data extraction was implemented in ArcGIS 10.2 (ESRI). Finally, the summary of environmental and richness data is shown in Table S2.

| Data analyses
We applied polynomial regressions (PRs) to estimate the relationship between species richness and elevation, guided by the corrected Akaike information criterion (AICc) value. We used generalized linear models (GLMs) to evaluate the elevational pattern for each ecological predictor, to assess the relationship between moss species richness and habitat diversity, and to predict the occurrence probability of each moss family along the elevational gradient by using presence/absence as the dependent variable. Acknowledged that these relationships may not be linear, especially that the elevational gradient of species richness can take many shapes but most often takes a hump-shaped pattern , we included polynomials of elevation up to the second degree for each GLM. We collected random samples from posterior distribution to estimate the 95% credible intervals for model parameters for the above models through the ARM package (Gelman & Su, 2016). The elevational trends of the predictors are shown in Figure 1. We fit multinomial models in a Bayesian framework using Markov chain Monte Carlo simulations in OpenBUGS through the package R2OpenBUGS (Sturtz et al., 2005), using elevation as the predictor and CT the outcome variable for the six most common moss families, respectively. Two Markov chains were simulated, each of length 10,000.
The burn-in was set to 1,000, and the chain was thinned by two to save work space and reduce autocorrelation. Convergence was assessed graphically and by the R-hat value (Brooks & Gelman, 1998).
We used an information-theoretic approach (Anderson et al., 2000;Chen et al., 2020;Stephens et al., 2005) to examine the relative roles of the MDE, the habitat complexity hypothesis, the energy hypothesis, and the environment hypothesis to moss species richness along the elevational gradient. Prior to analyses, all continuous variables (MDE, VT, CT, SR, AP, AT, and WS) were log-transformed to achieve normality and homoscedasticity. Then, we used GLM to build possible candidate models based on a priori hypotheses. Owing to the complication of habitat complexity hypothesis and environment hypothesis, we establish the models that included all the possible combinations of two habitat-related factors (VT and CT) for habitat complexity theory and three environmentrelated factors (AP, AT, and WS) for environment theory. A null model (richness ~ 1) was added for comparison. We calculated the variance inflation factor (VIF) of each variable in each model to assess collinearity. To reduce multicollinearity, only the models with VIFs < 10 were considered (Chen et al., 2020;Dormann et al., 2013).
We performed model averaging to evaluate the relative importance of each variable in shaping the elevational richness pattern (Galipaud et al., 2017). Last, we selected the best model through a forward stepwise selection algorithm. In the initial stage, we selected the model using a single variable with the minimum AICc. In the second stage, we then examined all two-variable models that included the variable chosen in the first step and chose the model with the minimum AICc. We then repeated the procedure for all three-variable models that included the two already selected, and so on, until AICc could not be further reduced. And this procedure was completed in the MuMIn package (Bartoń, 2015).
Thereafter, we used additive diversity partitioning to quantify the heterogeneity in biodiversity by comparing the diversity within (α) and among (β) sampling sites and by comparing the contributions made by elevation (β elevation ) and other factors (β replace ). In the additive approach, diversity can be explored across spatial scales (Gering & Crist, 2002), and γ-diversity (regional scale) is partitioned into the sum of the average diversity of sampling sites (α) and the heterogeneity among sampling sites (β). When a species is missing from a sampling site, one reason might be that the sampling site is bearing more geometric constraints of montane topography according to the MDE. So, we used additive diversity partitioning combined with species richness pattern predicted by the MDE and quadratic polynomial regression, respectively, to partition β into β elevation , which represents the average difference between α and the maximum diversity predicted by the MDE or quadratic polynomial regression (S max ) and β replace , the average number of missing species that are not explained by elevation. Because α, β, β elevation , β replace , and γ-diversity are measured using the same units, their relative importance can be quantified (Crist & Veech, 2006). We performed all analyses using R 3.5.2 (R Development Core Team, 2018).

| Species richness pattern
A total of 191 moss species, belonging to 73 genera under 26 families, were identified in 1,301 specimens at 73 sampling sites along the elevational gradient, in which four species (Drummondia sinensis, Sanionia uncinate, Tayloria indica, and Funaria hygromexrica) are endemic. Pottiaceae  (Table S3). Moss species in the Mt Tuofeng showed a hump-shaped richness pattern along the elevational gradient, with a distinct peak at 1,500 and 1,600 m a.s.l. (Figure 2). This result was con-  (Table 1). All four communities were recorded in five sampling sites with an elevation range between 1,338 and 1,678 m a.s.l. (Table S2), corresponding to the species richness peak.

| Relationship between species richness and explanatory factors
The information-theoretic statistics for the nine candidate models showed that MDE was suggested as the best model, which had an Akaike weight (W i ) of 0.98 and explained a largest proportion of variation for moss species richness pattern (R 2 = 0.245, p < 0.001; Table 2). Two alternative habitat-related models also provided a significant proportion of variation, one including only CT (R 2 = 0.146, p < 0.001) and the other including CT and VT (R 2 = 0.135, p < 0.01), but their ΔAICc exceeds two. The null model (richness ~ 1) had little support to the species richness pattern (ΔAICc = 19.33, W i = 0.00).
The best model for our dataset was the combination of MDE and CT, which reduced the AICc from 490.16 for MDE alone to 484.32 for MDE and CT, suggesting the MDE and the habitat complexity hypothesis were supported (Table 3).

| Additive partitioning of diversity
According to the additive diversity partitioning, α (12.79) explained only 6.7% of the variation in species richness, whereas β (178.21) explained about 93.3% of the variation in species richness (Figure 4).
We calculated the contribution of elevation toward the variation in species richness by using MDE prediction and quadratic polyno- We also compared β elevation and β replace through MDE prediction and quadratic polynomial regression, respectively. The proportion of β elevation to the total β ranged from 1.3% in the MDE prediction

| Species richness pattern of mosses on the Mt Tuofeng
The overall moss species richness pattern along the elevational gradient on the Mt Tuofeng is a hump-shaped pattern, peaking at midelevation between 1,500 and 1,600 m a.s.l., consistent with some studies in moss species (Grau et al., 2007;Wolf, 1993

TA B L E 2
Results of candidate models explaining variation for moss species richness pattern in the Mt Tuofeng other studies, moss species richness was found to have either no statistically significant trend Sun et al., 2013) or an increasing trend (Bruun et al., 2006) with altitude. The hump-shaped richness pattern is popular in studies from around the world, but there is no consistent explanation for this pattern. Our analysis showed that the MDE and the habitat complexity hypothesis concur to the elevational moss species richness. Model selection among alternative models showed that MDE is the primary driver for richness patterns, whereas little support was found for the energy and the environment.
It is not surprising that the energy hypothesis and the environment hypothesis are not supported in our data, as moss species are characterized by their poikilohydric condition and cold tolerance. The cuticle that seals the vascular plants body is often reduced or even lacking on the gametophyte of mosses, making mosses tolerant of desiccation and poikilohydric, which means that their water content is directly regulated by ambient humidity (Proctor et al., 2007). Moreover, a common feature among most mosses is their ability to grow at low temperature, and studies have showed that subglacial bryophytes following up to six centuries of ice entombment successfully regenerate (Cannone et al., 2017;La Farge et al., 2013;Roads et al., 2014). These ecophysiological features enable them to grow on rocks and tree trunks that are inhospitable for most vascular plants, thereby reducing the impact of energetic and environmental factors. The habitat complexity hypothesis was supported because habitat diversity plays an important role in maintaining biodiversity, and the removal of habitat types will obliterate species, especially habitat specialists (Sfenthourakis & Triantis, 2009). What's more, the CT pattern along the elevational gradient on the Mt Tuofeng is also hump-shaped (Figure 1b), implying a positive correlation between CT and moss species richness, which is proved by the GLM analysis (Figure 5b). A possible reason for the good support of MDE in our study is the high dispersal capacity of moss species (Patiño & Vanderpoorten, 2018), which makes a wide distributional range for most moss species, resulting in a high degree of overlap in the central area.
However, the drivers of elevational moss species richness pattern in our study are inconsistent with the studies of bryophyte diversity conducted by Raabe et al. (2010) and Spitale (2016) in European mountains where climate factors, such as temperature and SR, were the most important predictors. There are two reasons for this discrepancy. First, mosses and liverworts were sampled and analyzed in their studies. However, liverworts are usually more sensitive to drought than mosses (Oliver et al., 2000), so the drivers of elevational species richness pattern may differ between the two groups,  resulting in a disguised or biased pattern of mosses. Second, in their studies bryophytes were sampled from soil and wood, whereas we sampled mosses not only associated with soil and wood but also from dry rocks and stream water. Because humidity is higher on the forest floor than on tree trunks (Proctor & Tuba, 2002), bryophyte assemblages inhabiting deadwood and tree trunks are mostly subject to climatic variability (Spitale, 2016). Therefore, the sampling choice in our study could buffer the impact of climate.

| Three pieces of evidence support the MDE
In the tangled complexity of environmental and nonenvironmental factors affecting diversity gradients, new null models of the MDE helped to pare down the complexity, which is predicted where landmass boundaries such as oceans and mountaintops limit species ranges and the simple overlap of many, variously sized ranges create a peak in species richness at mid-elevation (Colwell & Hurtt, 1994;Colwell & Lees, 2000). In our dataset, we found three evidence conforming to the MDE. First, moss species richness peaked at mid-elevation ( Figure 2).
Second, within the 26 moss families, taxonomic groups that have a wider distribution usually peaked at mid-elevation, whereas those have a narrower distribution usually peaked somewhere away from the middle according to their occurrence probability along the elevational gradient ( Figure S1). Third, based on the edge effect or community overlap hypotheses, the greatest species richness exists in the ecotone areas of overlap between two distinct biological communities (Lomolino, 2001).
However, in our case, we found moss species richness had a negative relationship with the number of occupied VT (Figure 5a), implying the interior portion of the floral unit harbors the highest species richness, given the fact that the higher WS and lower humidity at ecotone areas are adverse to moss species survival (Liu et al., 2007). Therefore, our result, on the contrary, supports the MDE.

| How could elevation contribute to the heterogeneity
Elevation contributed toward explaining the heterogeneity, likely because of two facts. On one hand, the elevational gradient could provide heterogeneous environments (Hoorn et al., 2018;Rahbek, 1995), and different taxonomic groups survive at different elevation by selecting different physical conditions (Letten et al., 2013; Figure 6; Figure S1).
On the other hand, the elevational gradient may enable habitat segregation among moss species. Indeed, sympatric taxonomic species sharing similar resources should demonstrate some degree of niche overlap, leading to interspecific competition (Chesson, 2000;Dufour et al., 2015). In turn, to buffer competition and allow for coexistence, sympatric species may avoid each other in space and/or time and can generate differences in habitat selection (Holt, 1987;Milleret et al., 2018). Proportional use of CT along the elevational gradient varied among moss families (Figure 7), reflecting a certain degree of habitat segregation. This strategy of habitat use may reduce the effect of competition among different taxonomic groups. Although this phenomenon was found at the family level, it gives insights into the behavior of the 191 species in the system. And we speculate that habitat segregation among species may be more significant than among families. Both the two facts associated with elevation help to explain why species assembly varies among sampling sites, contributing the heterogeneity explained by elevation.  Figure 4). Moreover, elevation explained 2.2%, a weak predictive power, for the heterogeneity, and this may be due to four facts. First, according to the calculation formula (β elevation = S max − α), a lower predicted S max will further reduce the β elevation . Second, elevation or factors highly associated with elevation cannot adequately capture the high spatial heterogeneity of ecological and environmental variable characteristic of mountains. Third, mountain regions are home to aggregations of small-ranged species which could form centers of endemism. Fourth, moss species are supposed to have a high dispersal capacity; however, the large β replace in our study implies a high extinction rate when they colonize a new locality, as can be seen from Table S1 that many moss species appear discretely from sampling sites despite their wide distribution.

F I G U R E 5
To conclude, current species richness distribution pattern may bear the signatures of ecological and evolutionary effects, whereas evolutionary factors predominately shape the large heterogeneity through dispersal, extinction, and speciation processes. To further explore the extent to which each factor shapes the current pattern, we agree with

| Applicability to biodiversity conservation
The practical importance of these results for conservation is threefold. First, the positive relationship between species richness and CTs suggests habitat diversity is essential for sustaining species diversity, so conservation of habitat diversity is the key to maintain moss species diversity in the mountain. Second, the unimodal richness pattern we detected suggests that the highest moss species richness appears at mid-elevation; however, due to the large β-diversity and very small β elevation , conservation efforts should be paid to the whole elevational range rather than the mid-elevation only. Last, species of the moss families that have a very narrow distribution along the elevational gradient such as Drummondiaceae, Funariaceae, Splachnaceae, and Scorpidiaceae are likely to form small-ranged endemism ( Figure 6).
And special attention should be paid to preserve these irreplaceable species according to their distribution.

ACK N OWLED G M ENTS
This study was supported by the Natural Science Foundation of Hebei Province (D2019205019), the Department of Education

CO N FLI C T O F I NTE R E S T
The authors declare there are no competing interests.

DATA AVA I L A B I L I T Y S TAT E M E N T
All datasets used in this study are sourced from the literature which can be found in Supporting Information.