Evaluating the Scale Dependence of Biodiversity-productivity Relationships: A Case Study in Mountain Meadow Habitat

Ecological scale has been widely assumed to inuence various biodiversity-productivity relationships in ecological communities; however, its robustness has not been extensively studied. In this study, we tested the scale dependency of biodiversity-productivity relationships by evaluating their direct linkages while considering other confounders, and simultaneously incorporating functional traits and interspecic phylogenetic relationships. We surveyed sixty quadrats each with an area of 0.25 m 2 in three different meadows located along an elevational gradient in Yulong Mountain, China. We calculated different biodiversity parameters (richness, evenness, functionality, and phylogeny), and evaluated the chemical properties of soil from all quadrats, correlating their relationships with plant productivity at the local and regional scale. The direct and indirect relationships of biodiversity and productivity were evaluated using structural equation modeling. The biodiversity-productivity relationships were weak and inconsistent at the local scale, whereas some biodiversity metrics (richness, functional, and phylogenetic) showed either strong positive or negative relationships with productivity at the regional scale. However, a direct correlation between productivity and variables such as soil pH and community-weighted mean leaf carbon content was observed in the structural equation model reconstructed. Our study indicates that the scale dependency of biodiversity-productivity relationships in natural habitats may not be as strong as it may have been previously perceived, in case of taxonomic, functional, and phylogenetic diversity, respectively. Our study emphasizes the necessity to account for the confounding effects of abiotic factors when evaluating biodiversity-productivity relationship in natural habitats at regional or even worldwide scales.


Introduction
Plant biodiversity-productivity relationship as de ned in ecological research often represents the correlation between the number of plant species (i.e., taxonomic diversity) and the sum of their aboveground biomass production in a community (Whittaker and Niering 1975;Hector et al. 1999;Tilman et al. 2001). However, in context of the natural environment, the loss of species resulting in the potential loss of function of an ecosystem is highly debatable (Lavers and  The scale-dependency of taxonomic diversity-productivity relationships in ecological communities mainly result from distinct factors at the local and regional scales (Weiher 1999 . Plant taxonomic diversity-productivity relationships at local scales (among eld sampling plots in a relatively homogeneous environment) are expected to be positive when the number of species is low with a high probability of local species niche complementarity (whereby different species possesses different traits that utilize available resources differently) (Bond and Chase 2002). Whereas it is assumed to be negative when the number of species is high, with a high probability of species overlap in terms of resource utilization (Bond and Chase 2002). Regional scales when compared to local scales often contain distinct environmental gradients or ecological systems, which are often presumed to have stronger taxonomic diversity-productivity relationships due to their high environmental heterogeneitybased regional species diversity (Bond and Chase 2002). However, the strong correlation between taxonomic diversity and productivity at the regional scale might be a result of common responses to relevant environmental drivers or their internal linkages (Ma et al. 2010;Jing et al. 2015). A majority of studies conducted previously have focused on the correlation between taxonomic diversity and productivity without accounting for the effects of confounders such as abiotic factors (Ni et al. 2007; or have not explicitly evaluated the direct correlation between taxonomic diversity and productivity after controlling other relevant covariates (Maestre et al. 2012). Hence, a comprehensive understanding of the intensity of taxonomic diversity-productivity relationships at regional scales or local scales is lacking.
Compared to taxonomic diversity, functional diversity, which measures the mean or variance of single or multiple trait values of species in a community, might be more strongly related to productivity since functional traits take a direct part in the processes by which plants absorb, utilize, and transform essential resources (such as water and nutrients) ( Brun et al. 2019). Indeed, traits such as community-weighted mean (CWM), re ect the magnitude of the aggregation of trait values among species in a community, and higher CWM of functional traits may result in higher or lower productivity (Cadotte and Hillebrand 2017). Trait variance, which re ects functional complementarity (e.g., facilitation) or divergence (e.g., competition) among species in a community, may also have a positive or negative in uence on productivity (Cadotte and Hillebrand 2017).
Meanwhile, phylogenetic diversity, which is calculated using the phylogenetic relationship among species in a community is considered to be a parameter re ecting their functional similarity, and may have a closer relationship with productivity in comparison to taxonomic diversity (Cadotte et al. 2013). However, empirical studies to assess the relationship of functional and phylogenetic diversity with productivity at both local and regional scales in comparison with taxonomic diversity have not been reported till date.
In order to address the above-mentioned research questions, we surveyed 60 quadrats in three different meadow habitats, each measuring 0.5 m × 0.5 m, situated along an elevational gradient in the Yulong Mountains of Southwest China. For each species in the quadrats, we weighed their above-ground stem dry weight, while simultaneously recording a number of functional traits like maximum plant height (cm), leaf carbon, nitrogen, and phosphorus content (mg/g). A molecular phylogenetic tree for all the species in the quadrats was constructed using the DNA sequences of the coding genes rbcL and matK, respectively. For each quadrat, we also evaluated soil chemical properties such as pH, organic carbon (C, mg/g), total nitrogen (N, mg/g), and total phosphorus concentration (P, mg/g). In this study, we considered each meadow to be equivalent to the local scale, while all meadows taken together were considered to represent the regional scale. For each observed assemblage, we calculated the species richness, Shannon's evenness, and a series of functional and phylogenetic diversity metrics. Using multi-model inference and structural equation modeling, we aimed to test the following hypotheses: (1) plant taxonomic diversity-productivity relationship is stronger at the regional scale than at the local scale; however, the strong relationship at the regional scale might be driven by the responses of plant taxonomic diversity and productivity to the common abiotic covariates, and (2) functional and phylogenetic diversity have more robust relationships with productivity than taxonomic diversity over both local and regional scales.

Selection of study sites
Our study sites consisted of three different south-facing mountain meadows of regular topology along an elevational gradient (2,700,3,200 The mean annual precipitation is 935 mm, most of which occurs from June to October. The soil is typical upland red soil classi ed as wet iron bauxite (Gong 1999). There was an obvious species turnover in the plant communities of the meadows along the elevational gradient. The dominant genera were Isachne, Juncus at the lowest elevation, Ligularia, Agrostis at the middle elevation, and Agrostis at the highest elevation ). All the meadows have livestock-grazing histories (e.g., sheep, horses, and yaks).

Survey data
On each meadow, we randomly distributed twenty quadrats of dimensions of 0.5 × 0.5 m in August 2016.
The species composition within each quadrat was recorded, and the stems of each species of plant were cut at the ground level during the peak of the growing season. For many plants, it was di cult to distinguish individuals, so in those cases, we attempted to distinguish the individuals of the same species by separating their roots. All the stems at the ground level of each species were dried at 60 ℃ for 48 h, and then weighed to 0.1 mg. After cutting the stems at the ground level, we collected three different soil core samples in each quadrat using a cylindrical soil auger (5 cm inner diameter, 15 cm length). The soil cores from each quadrat were combined into a single composite sample, air-dried in the shade, and ltered using a 2-mm sieve for further stoichiometric analysis. Soil pH was measured from a 1:5 soil/water suspension. Soil organic carbon (C, mg/g), total nitrogen (N, mg/g), and total phosphorus (P, mg/g) were measured following the semimicro-Kjeldahl method, sulphuric acid-perchloric acid digestionbased method, and the Walkley-Black procedures, respectively (Institute of Soil Sciences 1978).

Plant traits and molecular phylogeny
Four plant traits that include maximum plant height (H max , cm), leaf carbon content (LC, mg/g), leaf nitrogen content (LN, mg/g), and leaf phosphorus content (LP, mg/g) were measured since these traits might re ect fundamental resource complementarity and interactions among co-occurring species Wright et al. 2004). We measured plant height for at most ve randomly selected samples of each species in each quadrat to account for the extent of the intraspeci c trait variability (de Bello et al. 2011). At least ve mature leaves per species in each quadrat were collected, dried to a constant weight at 60 ℃, and then weighted to 0.1 mg. The leaves per species were pooled together into a single composite sample across all quadrats. This was followed by assessment of leaf carbon, nitrogen, and phosphorus content using the same procedures used for analyzing the soil chemical properties.
A molecular phylogenetic tree was reconstructed for the plant species recorded in our study sites using the rbcL + matK regions of the chloroplast genome. The rbcL and matK sequences were aligned using MAFFT (Version 7.0, Katoh and Standley 2013) and then concatenated to form a supermatrix. For each gene, the top-ranked maximum-likelihood model of nucleotide substitution was selected using Akaike's information criterion, as implemented in the function "modelTest" in the phangorn library (Schliep 2011). A maximum-likelihood phylogeny was reconstructed using the BioNJ starting-tree using PhyML (Version 3.0, Guindon et al. 2010). One representative of early-diverging angiosperm lineage Amborrela trichopoda was chosen as the root of the phylogeny, which was then transformed to an ultrametric tree using the "chronopl" function with the parameter value of 1,000 in the APE library (Paradis et al. 2004).

Biodiversity assessment and determination of soil properties
For each assemblage, we calculated the species richness (S), which represents the number of species presented in an assemblage as the measure of the taxonomic diversity. We calculated the Shannon's evenness index (H′) with the abundance of each species (stem biomass production) in each assemblage using the function "diversity" in the vegan library (Oksanen et al. 2013). We also calculated a suite of single and multivariate functional diversity metrics based on plant traits, as well as phylogenetic diversity metrics using the maximum-likelihood phylogeny. The detailed descriptions for the measures of functional and phylogenetic diversity have been listed in Table S1. For functional diversity, we calculated the CWM for each single plant trait and multivariate functional diversity metrics such as functional dispersion (FDis), which sums weighted distances from a centroid in multivariate space (Laliberte and Legendre 2010). For the phylogenetic diversity, we calculated Faith's PD, which sums up the branch lengths of a phylogenetic tree linking all species of a particular assemblage (Faith 1992).

Data analysis
To determine the elevational patterns of plant productivity, the biodiversity metrics chosen for this study, and the soil chemical properties measured, we constructed a general or generalized linear model for each variable using the "lm" or "glm" function in the R stats library, respectively. For generalized linear models, the gamma distribution (e.g., productivity, Faith's PD, MPD, FRic, etc.) or the Poisson distribution (e.g., species richness) of the residuals was assumed to account for non-Gaussian error distributions and nonlinear relationships. The usage of a Gaussian distribution, a gamma distribution, or a Poisson distribution of model residuals was validated based on the normalized scores of standardized residual deviance (Q-Q plots). For each variable, we also constructed a corresponding intercept-only model as its "null" model that assumes no elevational pattern. Using both models for each variable, we calculated informationtheoretic ER, which explicitly evaluates the parameter bias-corrected likelihood of the null hypothesis against the alternative, as the sample size-adjusted Akaike's information criterion (AIC c , Burnham and To determine the relationships between plant productivity and the biodiversity metrics, as well as soil chemical properties measured at both local and regional scales, we constructed generalized linear models of plant productivity as a function of these standardized single predictors (mean/standard deviation) with the assumption of a gamma distribution of model residuals being similar for each meadow (local scale), as well as for all meadows as a whole (regional scale). For each predictor at each scale, a corresponding "null" model, which is the intercept-only model (where productivity ~ 1), was also constructed to calculate the ER. DE was calculated as a measure of the model's goodness-of-t.

Structural equation modeling
We constructed structural equation models (SEMs) to determine the strong correlation between plant productivity and individual biodiversity metrics at the regional scale driven by their responses to common relevant abiotic factors (such as elevation and soil chemical properties) or their internal linkages. Before model tting, plant productivity was log-transformed, and all strong predictor values (predictors with ERs > 3 at the regional scale) were standardized. The models were constructed and tted using the function "sem" in the library lavaan (Latent Variable Analysis) (Rosseel 2012). We started with a fully speci ed model whereby we assumed that: (1) plant productivity was dependent on a combination of the selected strong predictors, (2) plant biodiversity measures were directly linked with abiotic factors such as elevation and soil chemical properties, (3) functional and phylogenetic diversity measures relied on plant taxonomic diversity to test whether they mediate plant taxonomic diversity-productivity relationship, (4) soil properties were directly related to elevation if the elevation was pre-selected. Moreover, we assumed the correlations among functional and phylogenetic diversity metrics, as well as those existing among different soil chemical properties. The model constructed was then improved by eliminating the least nonsigni cant relationship, one at a time, until none remained. The statistical support for the nal model vs. the full model was assessed using AIC. The overall t of the nal model was evaluated using the comparative t index (CFI) and the root mean square error of approximate (RMSEA), with CFI ≥ 0.94 and RMSEA < 0.08 representing a good model t. Predictors' explanatory power was inferred from their respective R 2 value and standardized regression coe cient. All of the analyses were performed using R 3.5.3 (R Core Team 2015).

Relationships between plant productivity and biodiversity metrics and soil chemical properties at local and regional scales
To determine the scale-dependency of the relationships between the single predictors and productivity, we considered those showing strong relationships with productivity at the regional scale (i.e., the three meadows as a whole), including species richness (S), Faith's PD, MNTD, functional richness (FRic), CWM LC , soil pH, and total P for simplicity ( Fig. 2a-g). The most in uential predictor among these predictors was soil pH (ER = 10516, De = 29%), and followed by community weighted mean of leaf carbon content (CWM LC , ER = 469, De = 21%). However, no single predictor showed strong and consistent relationships with productivity at the local scale (i.e., each meadow) along the elevational gradient ( Fig. 2h-ab). The parameters which showed moderate relationships with plant productivity at the lowest and highest meadow include soil pH (ER = 3.31, De = 23%) and CWM LC (ER = 9.6, De = 30%), respectively.

Direct and indirect effects of biodiversity metrics on plant productivity at the regional scale
The nal SEM ( 2 = 21.837, p = 0.149, CFI = 0.988, RMSEA = 0.078, p = 0.269, Fig. 3) showed that both soil pH and CWM LC mediated the distribution of plant productivity along the elevational gradient, explaining > 38% of its variation. Species richness mediated the elevational patterns of MNTD (R 2 = 0.24), FRic (R 2 = 0.88), and PD (R 2 = 0.74). Interestingly, PD was additionally related to elevation although the link was quite weak. For total P, there was only a path from elevation.

Scale-dependency of explaining the variation in plant productivity
Our results showed that the taxonomic diversity-productivity relationship was only strong at the regional scale (ER = 212, De = 19%), which seems to support the hypothesis of scale dependency (Bond and Chase 2002). The taxonomic diversity-productivity relationship was the weakest at the mid-level elevation (De < 1%). However, functional evenness (FEve), functional divergence (FDiv), and functional dispersion (FDis), of which higher values may be re ective of higher e ciency of resource utilization and lower resource competition (Mason et al. 2005), also tended to be the highest at the middle elevation ( Fig. 1h-1j). This observation does not seem to corroborate with the prediction that taxonomic diversity-productivity should be positive at the local scale, when species in communities possess different traits that utilize available resources differently (i.e., complementary effect, Loreau and Hector 2001; Bond and Chase 2002). However, Cadotte and Hillebrand (2017) have indicated that stronger relationships between plant productivity and functional diversity metrics should exist when it is based on multiple traits than those based on single traits, when the complementarity effect drives the taxonomic diversity-productivity relationship. However, we found that the community weighted mean of single traits outperformed multidimensional trait measures in explaining the observed variation in plant productivity at all local scales (Table S2). Thus, our results may imply that selection effect whereby plant productivity is determined by dominant species (Loreau and Hector 2001), may be the main underlying mechanism in uencing the taxonomic diversity-productivity relationships at the local scale.
In our study we found that scale-dependency seemed to apply for functional/phylogenetic diversityproductivity-based relationship at least in case of strong regional predictors. This is in contrast with our expectation that the functional/phylogenetic diversity-productivity relationship would be more robust and more consistent than the taxonomic diversity-productivity relationship at both local and regional scales. A number of studies have reported the signi cance of functional/phylogenetic diversity over taxonomic diversity in explaining variation in plant productivity (Cadotte et al. 2009;Flynn et al. 2011;Liu et al. 2015Liu et al. , 2018. Our results showed that the top-ranked single models for all local scales that included either functional diversity or phylogenetic diversity generally agree with previous studies. However, our results also highlight the weak and inconsistent roles played by of functional/phylogenetic diversity in in uencing plant productivity at local scales. Our results may further suggest the dominant role of the selection effect for plant productivity at the local scale since the community weighted mean of single traits represented the top-ranked predictors for productivity in case of two of the three meadowlands studied. In addition, the relationships between plant productivity and soil chemical properties, which include soil pH and total P, appear to vary with spatial scale. It is not surprising that there was a strong relationship between plant productivity and soil pH at the regional scale since soil pH ranged from 4.28 to 6.14 along the elevational gradient (soil pH below 4 and above 8 is known to be detrimental to plant productivity, Sanchez et al. 2003) had the strongest elevational pattern (De = 68%) among the variables considered.
However, at the local scales, their relationship was signi cant in case of the lowest meadow considered in this study. One possible explanation is that the distribution of soil pH at the lowest meadow was more right-skewed compared to its approximate normal distributions at the middle and highest meadows (Fig.  S1). This may imply that plants at the lowest level prefer a neutral environment, whereas plants at the mid-level or at higher elevations have adopted to acidic soil conditions. Total P is the second most sensitive variable with respect to change in elevation, exhibiting a strong relationship with plant productivity at the regional scale, with no signi cant correlations at the local scale. Interestingly, there was a shift in the direction of the relationship between total P and productivity, and correlation between total P and soil pH transformed from negative to positive with change in elevation even though the concentration of total P was the highest at the highest elevation ( Fig. 2 and Fig. S2). This may imply that P is more important for plant production at higher elevations due to its functional role in promoting root growth and winter hardiness eliciting rapid gain of maturity, and stimulating tillering (Malhotra et al. 2018).

Strong biodiversity-productivity relationships at the regional scale do not guarantee presence of the internal linkages
The nal SEM results showed that there was no direct relationship between species richness and productivity, implying that their strong relationship at the regional scale might be driven by their responses to elevation. This also suggests that their relationship at the regional scale might not be as strong as it appeared after ruling out the in uence of elevation, which in turn might be consistent with their relationships at the local scale. To some extent, this reduces the strength of the scale dependency in taxonomic diversity-productivity relationships, and has important implications for other studies investigating taxonomic diversity-productivity relationships at the regional scale. Moreover, we found that the strong relationships between plant productivity and PD, MNTD, and FRic may be a result of their close connections with species richness. Interestingly, a direct link between PD and elevation was also discerned after accounting for the effect of species richness. This may be suggestive of the potential role of biogeographical processes (such as colonization, establishments, and extinction) for community phylogenetic diversity along the elevational gradient (Miao and Jianmeng 2015).
In this study, we found that both soil pH and community weighted mean of leaf carbon content (CWM LC ) mediated the distribution of plant productivity along the elevational gradient. Interestingly, pH of the soil and CWM LC showed the strongest elevational patterns (Fig. 1m and 1p). Furthermore, both were the only predictors having strong relationships with productivity at the local scale (Fig. <link rid=" g2">2</link>l and 2 aa). All these observations imply that a close connection between the local scale and regional scale may exist from the perspective of explaining the variation in plant productivity using only regional effective predictors (species richness, MNTD, PD, and FRic), as well as both local and regional effective predictors (soil pH and CWM LC ). Along with the insights provided by SEM in this study, it should also be noted that other variables not included in this study may potentially in uence the internal ecological relationships detected here, and hence, more in-depth studies are needed to get a more comprehensive understanding of these underlying ecological processes.

Conclusions
Scale-dependency has been widely reported in empirical studies predicting a stronger relationship between taxonomic diversity-productivity at a regional scale than the local scale (Craven et al. 2020). In this study, we found that the relationship between species richness and productivity was signi cant only along the elevation gradient, which is in line with previous studies (Ni et al. 2007;Li et al. 2019). However, the results of the SEM constructed in this study only supported the direct links between plant productivity and soil pH and CWM LC , which in uenced the productivity at both local and regional scales. This implies that the scale-dependency might not be as strong as previously reported. Our data supports the robustness of the impact of functional/phylogenetic diversity over taxonomic diversity in explaining the variation in plant productivity at both local and regional scales. Further evaluation of the direct linkages between local predictors and productivity using SEM at regional or even larger scales are pertinent for gaining in-depth understanding of these ecological processes.