Intraspeciﬁc Trait Variation Dilutes Deterministic Processes in Community Assembly of Arid Shrubs across Multiple Scales

: Trait-based approaches present a promising avenue for improving our understanding of species coexistence and community assembly, while intraspeciﬁc trait variation (ITV) across di ﬀ erent spatial scales is important in trait-based community assembly mechanisms, especially in extreme environments. In this study, we focused on the functional diversity and community assembly patterns of a desert community across di ﬀ erent spatial scales and investigated whether ITV plays a signiﬁcant role in community assembly processes in arid habitats. A 50 m × 50 m plot with di ﬀ erent small quadrats was established in a typical desert community at the transition zone between the Tengger Desert and Loess Plateau in China. A total of 14 traits were selected to assess the trait-based functional diversity and assembly processes in the community. We found that functional diversity showed di ﬀ erent patterns when considering ITV and related to di ﬀ erent types of traits (chemical traits or morphological traits) and some soil factors (pH and nitrate nitrogen). Plant communities in this study showed stochastic distribution patterns and similar functional diversity patterns based on functional trait approaches, regardless of spatial scales. Also, the e ﬀ ect of ITV on community assembly did not show more e ﬀ ect with increasing scales. These results indicated that ITV diluted deterministic processes in community assembly across scales in arid habitats.


Introduction
Community assembly is the process that shapes the structure and composition of plant communities, which describes the mechanism of diversity maintenance in a community [1,2]. Approaches based on functional traits have been widely used to demonstrate these assembly processes and diversity patterns [3][4][5][6][7][8]. However, the dominant trait-based community assembly or diversity studies have focused largely on trait differences among species [9], but there has recently been renewed interest in the role of intraspecific trait variation (ITV, [9][10][11]), which is extensive in nature and plays a significant role in trait-based community ecology [10][11][12]. A recent study by Siefert et al. (2015) estimated that an average of 25% of all trait variation in plant communities is intraspecific variation.

Study Area
The study area is situated in Jingtai County, Gansu Province, China (103 • 33 -104 • 43 E, 36 • 43 -37 • 38 N). The site is located at the transition zone among the Tengger Desert, Qilian Mountains, and the Loess Plateau [40]. This area has a temperate continental arid climate. The main geomorphological types are erosion-prone hills and piedmont pluvial and sloping alluvial plains. The mean annual precipitation is 183.5 mm, which mainly falls between June and September. The mean annual temperature is 8.1 • C. The mean annual total potential evapotranspiration is 3038 mm, and the frost-free period is 141 days [41]. The land surface is extremely dry, with severe desertification. Woody vegetation is dominated by small shrubs [42]. The dominant shrub species is Reaumuria soongorica, in the family Tamaricaceae.

Sampling Site Establishment
To examine the availability of sampling sites, we established three 50 m × 50 m sites across the Tengger Desert in Jingtai County (in the northern part of Jingtai County, the Tengger Desert stretches for about 16 13.80 N; Site RE-2 in Hongwawa at 103 • 58 06.6 E, 37 • 23 00.96 N. Each site is 7.8 km apart, we established three sites randomly within the target range of 100 m). We also obtained soil samples and documented the species information in preliminary experiments in the summer of 2016. In extremely dry deserts, habitats are largely homogeneous at the landscape scale (see Supplementary Materials Table S1). The pre-analysis proved that the three sites contain mostly redundant information because there were very similar habitat indicators (Supplementary Materials  Table S1). In addition, this study set up the biggest scale of 25 m × 25 m, the target species in this  Table S2, the target species had a relatively small individual canopy range in this desert habitat. In this study area, most of the woody species are undershrubs or subshrubs, with a height of only 5 cm and canopy diameter of less than 30 cm. Many more species are herbaceous, or small shrubs that appear only once. Thus, even within a 5 m × 5 m scale, the species richness or abundance that can be obtained is sufficient), so again, in the final analysis, one 50 m × 50 m site (104 • 08 15.55 E, 37 • 21 49.99 N, see Supplementary Materials Figure S1) was selected because this site has the greatest species richness, which can cover more regional species pool information (see Supplementary Materials Table S1, the richness of this site is 14 vs. 12 and 13 in RE-1 and RE-2). Although there was just one 50 m × 50 m plot, we treated it as a sampling site area or a frame of reference rather than a repeatable plot. We divided the plot into non-overlapping square quadrats of 5 m × 5 m, 10 m × 10 m, 15 m × 15 m, 20 m × 20 m, or 25 m × 25 m scales (25,100,225,400, and 625 m 2 , respectively; Figure 1). The number of each scale plot was 100, 25, 5, 4, and 4, respectively. Community investigation was confined to the shrub and herbaceous layer. Each species in 50 m × 50 m was identified and documented with its abundance, coverage, plant height (H), and canopy range (CR, calculated by length and width). Also, the abundance and coverage of species in 5 m × 5 m, 10 m × 10 m, 15 m ×15 m, 20 m ×20 m, and 25 m × 25 m were documented. In total seven shrub species with 2412 individuals and seven herbaceous species with 3839 individuals were documented (See Supplementary Materials Table S2). Supplementary Materials Table S1). The pre-analysis proved that the three sites contain mostly redundant information because there were very similar habitat indicators (Supplementary Materials  Table S1). In addition, this study set up the biggest scale of 25 m × 25 m, the target species in this  research area (50 m × 50 m) are enough (see Supplementary Materials Table S2, the target species had a relatively small individual canopy range in this desert habitat. In this study area, most of the woody species are undershrubs or subshrubs, with a height of only 5 cm and canopy diameter of less than 30 cm. Many more species are herbaceous, or small shrubs that appear only once. Thus, even within a 5 m × 5 m scale, the species richness or abundance that can be obtained is sufficient), so again, in the final analysis, one 50 m × 50 m site (104°08′15.55″ E, 37°21′49.99″ N, see Supplementary Materials Figure S1) was selected because this site has the greatest species richness, which can cover more regional species pool information (see Supplementary Materials Table S1, the richness of this site is 14 vs. 12 and 13 in RE-1 and RE-2). Although there was just one 50 m × 50 m plot, we treated it as a sampling site area or a frame of reference rather than a repeatable plot. We divided the plot into nonoverlapping square quadrats of 5 m × 5 m, 10 Table S2).

Soil Sample Characteristics Measurement
We collected soil samples (three repetitions from different aspects of each hole) from 0-20 cm below the ground from the center of each 5 m × 5 m sub-plot. After sampling, we analyzed the soil water content (SWC), pH, and contents of ammonium nitrogen (AN), nitrate nitrogen (NN), and rapidly available phosphorus (RAP) based on the average of three soil sample repetitions to determine the impact of soil character on the community. For 10 m × 10 m, 15 m × 15 m, 20 m × 20 m, 25 m × 25 m scales, these soil characters were averaged by 5 m × 5 m sub-plots at four corners of each plot. These stoichiometric factors of soil are measured by a chemical analyzer (Cleverchem 200+, EuroVector Inc., Pavia, Italy) and a spectrophotometric method. We also recorded the altitude, longitude, and latitude of the soil sampling site, which were all obtained by a GPS recorder (ProMark 100, Ash Tech, Inc., Winthrop, MA, USA).

Plant Sampling and Trait Measurement
Research on functional trait needs to maximize the dimension of traits but minimize the number of traits. Thus, traits belonging to the leaf, stem, and root should be contained (Laughlin 2014). In this study, both morphological and stoichiometric traits were considered. We recorded the plant height (H) and coverage (CR) of every woody individual occurring in the plots. For each woody species occurring in the 5 m × 5 m sub-plots, we collected no less than 20 fully mature leaves from as many individuals as possible. Stems and roots of woody species were collected from each sampled individual. For each herbaceous species, we collected 5-10 whole plants and then divided them into leaf, stem, and root sections.  Figure 1) was resampled according to the above standard. Every leaf we collected was scanned by Epson Perfection (V-700 Photo, Epson Inc., Nagano, Japan). The leaf area (LA) and leaf length (LL) were measured by Motic Images Plus 2.0 (Motic Digital Microscopy Co., Xiamen, China). In addition, we measured the specific leaf area (SLA), leaf dry matter content (LDMC), stem dry matter content (SDMC), and root dry matter content (RDMC) following standard methods [43]. We also measured the stoichiometric traits of plant samples. The leaf nitrogen content (LNC), leaf carbon content (LCC), stem nitrogen content (SNC), stem carbon content (SCC), root nitrogen content (RNC), and root carbon content (RCC) were measured by the elemental analyzer (EA3000, EuroVector, Inc., Pavia, Italy), whereas the leaf phosphorus content (LPC) was measured by CleverChem 2000 (Dechem-Tech, Inc., Hamburg, Germany) using standard methods [44]. In addition, we established a "multiple-trait" [45,46], which may reflect the response of the community traits to the environments directly as it was a unidimensional value. To avoid collinear redundancy, each above trait underwent principal component analysis (PCA, see Supplementary Materials Table S2). Finally, we selected LA, LDMC, RDMC, LNC, SNC, RCC, and LPC to calculate the multiple-trait value of one community by multidimensional hypervolume algorithm based on Euclidean distance among these traits [47].

Functional Diversity
We calculated Mason's functional α-diversity (F α ) [48] and Mason's functional β-diversity (F β ) [49]. The functional diversity of each scale of communities with and without considering for ITV (method with/without ITV, see [50][51][52]). Functional diversity was calculated by using the "FD" package in R 3.6.1. The algorithms of F α and F β , are as follows: Diversity 2020, 12, 447

of 14
Without considering ITV, P i is the relative abundance of species i, x i is the functional trait value of species i, x is the mean functional trait value of all species, x k is the community weighted mean trait value of community k, and x region is the mean functional trait value of the whole study region. Considering ITV, P i = 1, x i is the functional trait value of individual i, x is the mean functional trait value of all individuals, x k is the mean functional trait value of individuals in community k, and x region is the mean functional trait value of the whole study region. For plots at each scale, we calculated the functional diversities both with and without ITV and analyzed the significance of difference values between with/without ITV by paired t-test. The t-test was run on the STATISTICA 6.0 Program. All scatter diagrams were constructed with Origin Pro 8.0 software.

Functional Trait-Based Community Assembly Pattern
At each sampled spatial scale, we calculated the standardized effect size (SES) which represents the deviation of traits with/without considering ITV [52]. The SES described the level of trait convergence or divergence compared to the null model. We calculated these values applied in the R-package "Picante" with 999 repetitions [31,53]. Note that the SES values of traits were weighted by each species' relative abundance. A negative SES (especially SES < −1.96) value indicates that the functional diversity is lower than the expected random distribution, while a positive SES (especially SES > 1.96) value indicates that the functional diversity is higher than the expected random distribution [54]. For plots at each scale, we calculated the SES both with and without ITV and analyzed the significance of the difference between values with/without ITV by paired t-test. The t-test was run at STATISTICA 6.0 Program (Statsoft Inc., Las Vegas, NV, USA). All scatter diagrams were constructed with Origin Pro 8.0 software (OriginLab Co., Studio City, CA, USA).

Effect of Environmental Conditions on Functional Diversity
To reveal whether extreme environmental conditions affect the diversity and community assembly patterns, we used multiple regression considering alpha and beta diversity as response variables and the first axis of the PCA including soil physical and chemical variables (e.g., AN, NN, pH, RAP, and SWC) as factors. A PCA of the five soil traits (including AN, NN, pH, RAP, and SWC) was performed to represent the complex soil factors at five spatial scales. As response variables, we used the functional alpha and beta diversity indices at each sampled spatial scale. A PCA was conducted with SPSS 19.0 software (IBM, Armonk, NY, USA).

Functional Diversity and Community Assembly Patterns
When considering ITV, there is a higher functional alpha diversity in multiple-traits and most single-traits (such as LCC, LDMC, LPC, LNC, SDMC, SCC, SNC, RCC, RNC, RDMC, and CR) than that in the analysis without considering ITV ( Figure 2). Moreover, our results from analyses considering ITV showed lower functional alpha diversity for some single-traits (such as LL, SLA, and H) than that in analyses without ITV ( Figure 2).
We found that considering ITV would result in higher functional beta diversity both in multiple-trait and some single-trait analyses (such as LNC, LPC, LCC, SDMC, SNC, SCC, RCC, RNC, and CR) than that in the analysis without ITV (Figure 3). On the other hand, our results from analyses considering ITV showed lower functional beta diversity for a few single-trait analyses (such as SLA, LL, and H) than that in the analysis without ITV. In addition, some single-trait analyses (such as LDMC and RDMC) showed no significant difference for functional beta diversity with considering ITV to those without considering ITV (Figure 3). We found that considering ITV would result in higher functional beta diversity both in multipletrait and some single-trait analyses (such as LNC, LPC, LCC, SDMC, SNC, SCC, RCC, RNC, and CR) than that in the analysis without ITV (Figure 3). On the other hand, our results from analyses considering ITV showed lower functional beta diversity for a few single-trait analyses (such as SLA, LL, and H) than that in the analysis without ITV. In addition, some single-trait analyses (such as LDMC and RDMC) showed no significant difference for functional beta diversity with considering ITV to those without considering ITV (Figure 3). The SES of traits with ITV/without ITV were significantly different from each other. We also found that considering ITV could produce a higher SES for both multiple traits and some single traits (such as LCC, LL, LPC, LDMC, RCC, and RDMC) than those without ITV. However, our results also indicated that considering ITV could produce a lower SES for a few of the single traits (such as LNC, RNC, SCC, SNC, and CR) than those without ITV. In addition, some single traits (SLA, SDMC, and H) showed no significant trend for SES when comparing analyses with ITV to those without ITV (Figure 4). In addition, when considering ITV, SCC showed a significant divergent assembly pattern (SES > 1.96), while LPC showed a significant convergent assembly pattern (SES < −1.96).

The Effects of the Soil Environmental Factors on Functional Diversity
The effect of environmental conditions on the functional alpha diversity was analyzed for multiple soil traits (PCA1, see Supplementary Materials Tables S3 and S4) and each single soil trait (such as AN, NN, pH, RAP, and SWC) at each spatial scale (Supplementary Materials Figure S2). There was a negative correlation between PCA1 and the functional alpha diversity with/without ITV at the Diversity 2020, 12, 447 7 of 14 10 × 10 m 2 spatial scale, as was observed with pH. Additionally, there was a negative correlation between NN and functional alpha diversity with ITV at the 10 × 10 m 2 spatial scale only. At the other spatial scales, we found no significant correlation between environmental heterogeneity and functional alpha diversity (Supplementary Materials Figure S2). The SES of traits with ITV/without ITV were significantly different from each other. We also found that considering ITV could produce a higher SES for both multiple traits and some single traits (such as LCC, LL, LPC, LDMC, RCC, and RDMC) than those without ITV. However, our results also indicated that considering ITV could produce a lower SES for a few of the single traits (such as LNC, RNC, SCC, SNC, and CR) than those without ITV. In addition, some single traits (SLA, SDMC, and H) showed no significant trend for SES when comparing analyses with ITV to those without ITV (Figure 4). In addition, when considering ITV, SCC showed a significant divergent assembly pattern (SES > 1.96), while LPC showed a significant convergent assembly pattern (SES < −1.96). We found weak correlations between functional beta diversity and some soil factors with or without ITV at the 5 m × 5 m spatial scale (e.g., RAP, pH and AN; Supplementary Materials Figure S3) and 10 × 10 m 2 spatial scale (e.g., PCA1, NN, AN, SWC, and RPC; Supplementary Materials Figure S4). However, for other scales, we found no significant correlations between environmental heterogeneity and functional beta diversity.

Functional Diversity and Community Assembly Patterns across Scales
The effects of spatial scale and ITV on the functional diversity within microsites (alpha diversity) were analyzed for multiple traits and every single trait at different spatial scales with and without ITV (Figure 2). Functional alpha diversity had no significant differences across each spatial scale with or without ITV. The relationship between spatial scales and functional alpha diversity significantly differed between that with ITV and that without ITV for most functional traits (Figure 2).

The Effects of the Soil Environmental Factors on Functional Diversity
The effect of environmental conditions on the functional alpha diversity was analyzed for multiple soil traits (PCA1, see Supplementary Materials Tables S3 and S4) and each single soil trait (such as AN, NN, pH, RAP, and SWC) at each spatial scale (Supplementary Materials Figure S2). There was a negative correlation between PCA1 and the functional alpha diversity with/without ITV at the 10 × 10 m 2 spatial scale, as was observed with pH. Additionally, there was a negative correlation between NN and functional alpha diversity with ITV at the 10 × 10 m 2 spatial scale only. At the other Our results showed that the functional beta diversity did not significantly differ across each spatial scale with or without ITV (Figure 3). The observed functional beta diversity, the relationship between Diversity 2020, 12, 447 9 of 14 spatial scale and functional beta diversity, significantly differed between analyses with ITV and those without ITV for most functional traits (Figure 3).
The SES of functional traits did not significantly differ across each spatial scale with or without ITV (Figure 4). Almost all SESs of traits showed random assembly patterns across scales, regardless of considering ITV or not. SES of SLA showed different patterns at 15 m × 15 m, 20 m × 20 m scales when comparing to other scales. SES of SDMC showed different patterns at 10 m × 10 m, 15 m × 15 m scales when comparing to other scales.

Discussion
As expected, our study emphasizes the importance of considering ITV when studying functional diversity patterns in an arid ecosystem with species-poor communities. In addition, our results showed that ITV could influence the observed functional diversity patterns [55][56][57]. As well as at small scales at which plant species manifest with competitive and facilitative interactions through phenotypic plasticity.
Some previous studies considered that the importance with/without ITV was highly dependent on the spatial scale [23,26]. A previous study has found in tropical forests that both habitat associations and niche differentiation shape species co-occurrence patterns from small to intermediate spatial scales [19]; Münkemüller et al. [58] suggested that scale choice is decisive in revealing signals in diversity patterns in European Alps. Nevertheless, these patterns show inconsistent results from our study. Our results indicated that ITV significantly affects functional diversity at any spatial scale, and considering ITV for every single trait does not lead to consistent results within plant communities. Functional alpha diversity based on leaf "chemical" traits tended to have higher alpha diversity when considering ITV in plant communities, which agrees with previous studies that found a high ITV in leaf chemical traits [21,26,59]. In contrast, morphological traits with ITV (such as LL, SLA, and H) showed lower functional alpha diversity than those without ITV (Figure 2), which may be due to high ITV. The pattern of functional beta diversity is similar to that of functional alpha diversity; beta diversity was higher based on chemical traits but lower based on morphological traits with ITV, suggesting a consistent trend for functional alpha and beta diversity beyond scales. In addition, among spatial scales, the functional alpha and beta diversities hardly shifted, suggesting that the traits that we selected were conservative; otherwise the habitat conditions were too homogeneous across scales based on our sampling design, which was confirmed by our pre-analysis.
In addition, environmental factors could also influence functional diversity, for example, Shiono et al. [60] reported that species richness decreases with increasing climate harshness in temperature and precipitation seasonality zones. The extremely dry conditions in our study area led to distinctive patterns of functional diversity. In our study, we measured five soil factors as direct environmental factors. We found that functional alpha diversity decreased as the variation in soil PC1 scores and NN increased, whereas alpha diversity increased along with pH (Supplementary Materials Figure S2). However, functional beta diversity changed weakly along several soil variation gradients (Supplementary Materials Figures S3 and S4). For other soil factors, there were no significant trends with functional beta diversity, indicating very weak relationships between functional beta diversity and soil factors resulting from the homogeneous conditions in our study.
In most previous studies, researchers generally paid great attention to the community assembly mechanism in tropical forests, which have high species richness, few restrictive factors, and rare niche differentiation [61,62]. The pattern of community assembly influences the diversity pattern. Tirado et al. [63] found through investigating the neighbor effect that positive interactions among species prevail in communities where positive species associations dominate; the appearance of benefactor species in patches increases species richness compared with that in the surrounding inter-shrub spaces in different semi-arid habitats in southeast Spain. Chalmandrier et al. [21] found that, at large spatial scales, the importance of the imprint of environmental filtering on functional diversity increases with increasing stress, while the effect of intraspecific variability on assembly rules is noticeable only at a small spatial scale in grassland communities [21]. In this study, a null model approach was used to test for trait-based assembly processes of shrub communities in an arid environment. When considering ITV, dissimilar trait-based assembly patterns would be shown based on different traits, in order to adapt to some environmental factors in an extreme arid condition. Here, LCC and RCC were convergent considering ITV, while SCC was divergent. Without considering ITV, filtering processes based on traits related to carbon (except SCC) were observed, while divergence patterns based on all traits related to nitrogen were observed, which may be caused by species interaction or even environmental filtering [64][65][66]. These results suggested that the carbon content was relatively conservative between different species, whereas the nitrogen content showed great variation since nitrogen participates in many plant metabolic processes. The SES based on the dry matter content (except SDMC) also showed filtering processes, indicating that desert plants have a tendency to adapt to the environment under water stress (dry matter content and moisture content are weighted against each other). The patterns of SCC and SDMC were obviously different from those in leaves and roots, which indicated that the stem traits of desert plants may be different from traits of leaves and roots, which is worth further study. Without considering ITV, the SES of multiple traits showed filtering patterns at different spatial scales. This result indicated that the extreme environment in the desert would result in a more consistent combination of traits, that is, environmental filtering may dominate the community assembly mechanism in this study area.
Besides, we found no significant change tendency in relationships of diversity and assembly pattern when considering ITV and without considering ITV across different scales. It indicated that the effect of ITV which influencing functional diversity and community assembly would not change with increasing scales.
For a community assembly study, there may be different results with or without considering ITV, suggesting that ITV is important for studying assembly patterns [10,11,24,26]. On the other hand, since more data enter the functional structure model when considering ITV, considering ITV seems more realistic than not considering it. Interestingly, in our study, almost all traits (including multiple-traits) exhibited a stochastic distribution when considering ITV. In other words, when considering ITV, community functional trait distributions may shift towards the opposite direction regardless of whether these traits have convergent or divergent patterns comparing to those without considering ITV. This phenomenon may lead to a consistent pattern of SES of these traits. Our results may draw a more "stochastic" pattern when considering ITV. Therefore, the neutral-like assembly pattern may be diluted by ITV, rather than by ecological processes. In general, in extreme environments, we used to expect the assembly pattern to be deterministic, but the current results suggest that the assembly pattern becomes more stochastic when ITV is considered. These results suggest that environmental filtering screens the weighted mean of community traits, but not individual trait values. Under extreme conditions, individuals still need to generate enough intraspecific variation to better occupy a limited niche. Therefore, when ITV is considered, traits that have been screened in the community may be masked, and the variation of traits displayed by individuals may lead to stochastic assembly patterns. Based on our findings, it is necessary to consider ITV when studying community assembly in arid regions in the future. Considering ITV can provide more information and a new perspective for studying the succession and development of communities in arid regions.

Conclusions
Our study in an extremely arid environment demonstrates the importance of ITV in diversity and community assembly analyses based on functional traits. Considering ITV greatly impacts functional diversity, and this effect is related to different types of traits (chemical traits or morphological traits). Specifically, our results indicate that the ITV effect does not change with increasing scales in community assembly. In addition, an extreme environment shapes specific combinations of traits. Our results also show that the community assembly process may be neutral at some levels and requires further study. When considering ITV in our study, the community assembly process may have been diluted, suggesting that ITV is important for studying assembly mechanisms in species-poor communities.
Supplementary Materials: The following are available online at http://www.mdpi.com/1424-2818/12/12/447/s1, Figure S1: Landscape vegetation distribution, Figure S2: Influence of environmental heterogeneity on functional alpha diversity, Figure S3: Influence of environmental heterogeneity on functional beta diversity, Figure S4: Influence of environmental heterogeneity on functional beta diversity, Table S1: Soil factors and species richness of each plot, Table S2: The species of this study, Table S3: PCA loadings of each functional trait on PCA 1 axis and PCA 2 axis, Table S4