Covariation in Plant Functional Traits and Soil Fertility within Two Species-Rich Forests

The distribution of plant species along environmental gradients is expected to be predictable based on organismal function. Plant functional trait research has shown that trait values generally vary predictably along broad-scale climatic and soil gradients. This work has also demonstrated that at any one point along these gradients there is a large amount of interspecific trait variation. The present research proposes that this variation may be explained by the local-scale sorting of traits along soil fertility and acidity axes. Specifically, we predicted that trait values associated with high resource acquisition and growth rates would be found on soils that are more fertile and less acidic. We tested the expected relationships at the species-level and quadrat-level (20×20 m) using two large forest plots in Panama and China that contain over 450 species combined. Predicted relationships between leaf area and wood density and soil fertility were supported in some instances, but the majority of the predicted relationships were rejected. Alternative resource axes, such as light gradients, therefore likely play a larger role in determining the interspecific variability in plant functional traits in the two forests studied.


Introduction
The distribution of species and communities along environmental gradients is a central focus in ecology. The distribution of species is expected to be determined by the distribution of resources. The functional strategy of a species will dictate its resource use and therefore its location along a resource axis or resource axes. Thus function should vary predictably along these gradients. This has lead to a tradition in plant ecology of predicting and analyzing the geographic distribution of functional strategies [1,2].
The relationship between plant function traits and environmental gradients has been quantified for a number of plant traits using large-scale datasets. Evidence from these broad-scale functional trait analyses suggest that the mean functional trait value of an assemblage changes predictably along environmental gradients. For example, leaf and wood traits, seed mass and maximum height have been shown to vary predictably with mean annual temperature [3][4][5][6][7][8][9][10]. Additional studies have also examined the relationship between leaf and wood traits with soil nutrient levels. Leaf economic traits related to resource acquisition such as specific leaf area, leaf nitrogen content and leaf phosphorus are positively correlated with soil nutrient content and these relationships were stronger than those with climatic gradients [11]. Wood density, which is negatively correlated with volumetric growth rates, is negatively correlated with nitrogen and phospho-rus levels across the Amazon Basin [12]. A running theme in many of these papers is there is a trade-off between the structural allocation and demographic rates based on the resource availability. Specifically, species that favor high resource environments should have higher growth and mortality rates where biomass is allocated to producing a large amount of small seeds that germinate quickly, structurally cheap leaves that have high specific leaf areas but photosynthesize at a high rate, and structurally cheap wood that permits rapid volumetric growth into the canopy. In contrast, species that favor low resource environments should be characterized by ecological strategies that increase structural investment at the cost of decreased resource acquisition and demographic rates. While many of the above studies have supported the expected relationships between environmental gradients and plant traits across broad gradients, this work has also demonstrated that a tremendous level of interspecific variation occurs within locations along the gradient [13].
The large inter-specific trait variation within sites in global datasets could be the result of trait -environmental gradient relationships on local scales and how different ecological strategies related to resource acquisition and demographic rates sort out along important resource axes. For example, given the previous research showing strong and consistent relationships between plant traits and soil nutrients on global scales, it is expected that local scale plant trait distributions should also vary predictably along local scale soil nutrient gradients. In particular, we predict that individual species with plant traits associated with high rates of resource acquisition and growth such as high values of specific leaf area, maximum height and leaf area and low values of seed mass and wood density are predicted to occur on soils with high nutrient content. Conversely, species with low values of specific leaf area, maximum height and leaf area and high values of seed mass and wood density are expected to be located in soils with low nutrient levels.
Here, we integrate tree distribution and soil nutrient data with five plant functional traits -specific leaf area, maximum height, leaf area, seed mass and wood density to test the predicted relationships among local-scale gradients in soil nutrient levels. In particular, we quantify: (1) the correlation between species mean trait values and their mean position on soil nutrient gradients and (2) the correlation between the mean trait value in 20620 m quadrats and the soil nutrient level in that quadrat. The analyses are performed separately in two forest inventory plots. The two forest plots were chosen for two important reasons. First, they share similar forest inventory, trait collection and soil nutrient mapping protocols making a comparative study feasible. Second, the forests are vastly different in their topographic heterogeneity thereby allowing us to determine whether the degree of local habitat heterogeneity influences the strength of trait-soil relationships. We first test the above predictions using species-level data and then ask whether the species-level relationships scale-up to the quadrat-level where the mean trait value within a quadrat can be predicted based on the soil nutrient levels in that quadrat.

Research Sites
The datasets used in this study were compiled from two permanent large forest dynamics plots in tropical and subtropical forests. The Barro Colorado Island (BCI) 50-ha forest dynamics plot (9u109N, 79u519W) is located on well-weathered kaolinitic Oxisols on Barro Colorado Island, Panama (Fig. 1), and it is characterized as a lowland semideciduous moist forest. In the 10006500 m rectangular area, the plot spans an altitudinal range of 120 to 160 m and the slope ranges from 7u to 20u. Daily maximum and minimum temperatures average 30.8uC and 23.4uC, respectively. Annual rainfall averages 2600 mm, with just 10% of the annual total falling during a 4-mo dry season. The BCI plot was first censused in 1981/82 [14]. All trees with a diameter at breast height (dbh) $1 cm were measured, identified and mapped. A second census was performed in 1985 and censuses have been repeated every 5 years thereafter. Here we use the 2005 census, which includes 208,387 individual trees belonging to 299 species.
The Gutianshan (GTS) 24-ha permanent plot (29u159N, 118u079E) is located in the old-growth forest of Gutianshan National Nature Reserve, Kaihua County, Zhejiang Province, Southeast China (Fig. 1), and it is characterized as a subtropical evergreen broad-leaved forest. The GTS forest plot contains approximately 140,000 individual trees (dbh$1 cm), representing 49 families, 103 genera and 159 species in the plot. It was established in the summer of 2005, following the same protocol as for BCI [15,16]. In the 6006400 m rectangular area, the plot spans an altitudinal range of 446.3 to 714.9 m and the slope ranges from 13u to 62u. The mean annual temperature in the Gutianshan Reserve is 15.3uC. The hottest month is July (mean temperature of 27.91uC), and the coldest is January (mean temperature of 4.31uC). Annual precipitation averages 1963.7 mm, with a dry period between October and February.
The major soils can be classified into four types: red soil, redyellow soil, yellow-red soil and marsh soil [15].

Plant Functional Traits
We measured leaf area (LA), specific leaf area (SLA), wood density (WD), seed mass (SM) and maximum height (H max ) for species at both sites. The trait collection protocols for BCI are described in Wright et al. [17], and the GTS collection protocols followed Cornelissen et al. [18] with the exception of WD which followed the protocols of Wright et al. [17]. Below we briefly describe the collection methods and sample sizes for the GTS plot.
Leaf area and SLA were measured using at least ten mature leaves collected from the tallest portion of the canopies of 5-10 of the largest individuals of each species. The SLA was calculated as mean of fresh leaf area divided by the leaf dry mass without the petiole. The LA was measured as the mean leaf surface area without petioles for each species. The SM was calculated by collecting 30 to 200 mature, fresh seeds from more than five individual trees of each species in or near the plot. We removed appendages and oven dried seeds for 48 h at 80uC. The SM value is the mean value over all seeds of each species. The WD was calculated by collecting wood samples from 5 to 10 individuals for each species in the area surrounding the plot using methods described in Wright et al. [17]. The H max values for GTS were estimated using values reported in the Flora of China [19] and the Flora Reipublicae Popularis Sinicae [20].

Phylogenetic Trees
Two phylogenetic trees were utilized in this study. Specifically, we utilized a phylogenetic tree from Kress et al. [21] for the BCI plot. This phylogeny was constructed using a DNA supermatrix composed of three sequence regions -rbcL, matK, and trnH-psbA. The supermatrix and the software RAxML [22] were used to construct a maximum likelihood phylogenetic tree. We constructed a phylogenetic tree for the GTS forest plot species following the same methodology as Kress et al. [21].

Soil Fertility
Soil samples were collected and analyzed following the protocols established by the Center for Tropical Forest Science (CTFS) in both plots (http://ctfs.si.edu/datasets/bci/soilmaps/BCIsoil.html) [23]. However, the sampling design and intensity differed. At BCI, the plot was divided into a 50650 m grid, grid intersections were basal collection points, and additional collection points were marked at 2 m, 8 m or 20 m along a random compass direction from each basal point. Thus, 300 points were sampled in the 50ha plot. At GTS, the grid was 30630 m, the additional collection points were 2, 5 and 15 m along a random compass direction from each basal point, and a total of 892 samples were collected inside the 24-ha plot.
John et al. [23] describe the methods used to process BCI soil samples. At GTS, a 300-400 g topsoil sample was taken from 0-10 cm depth and air-dried. The soil was then sieved with a 2 mm mesh screen. The sieved soil was used to extract available cations. 50 g of the 2 mm-filtered soil was filtered again with 0.15 mm mesh screen for analyses of total C, N, and P. Additional samples from a depth of 15 cm were taken using two polyethene pipes with a diameter of 5 cm. One of these samples was used for extracting NH 4 + and NO 3 2 (using 2.0 M KCl on 2 g soil) and measuring gravimetric moisture content and pH value (soil: water was 1:5 Additional detailed information regarding soil data collection can be found in John et al. [23] for BCI and Liwen et al. [24] for GTS. For each forest plot, we used a principal components analysis (PCA) to extract orthogonal axes of soil fertility and acidity from the 13 measured soil nutrients and to reduce information redundancy. We used the significant PCA axes to characterize soil fertility and acidity for all subsequent analyses.

Statistical Analyses
Our datasets included mean trait values (T) for each species S (T S ), soil fertility (F) for each 20620 m quadrat Q (F Q ) and the number of individuals of each species in each quadrat (N SQ ). We used these measured values to calculate mean soil fertilities for each species (F S ) and mean trait values for each quadrat (T Q ) as follows: We first performed a species-level analysis to test our predictions by calculating a Pearson correlation of trait values and mean soil properties for species calculated by weighting the PCA scores of each 20620 m quadrat by the number of individuals of species in that same quadrat (eqn. 1) The LA and SM values for species in both forest plots were log transformed to satisfy the normality assumption. Next we used phylogenetically independent contrasts (PICs) [25][26][27] to evaluate relationships between measured values of T S and calculated values of F S . This second analysis was used to factor out the bias of phylogenetic non-independence and to evaluate the hypothesis that evolutionary changes in trait values were associated with the spatial distribution of species with respect to soil fertility. PIC regressions were forced through the origin [28] and significance was evaluated after removing extreme outliers (absolute value of studentized residual.5) and contrasts with undue leverage (leverage.0.2).
We performed a third correlation analysis to evaluate the relationship between calculated T Q (eqn. 2) and measured/ estimated F Q . The LA and SM values from the BCI plot were log transformed to satisfy the normality assumption. This quadratlevel analysis was used to test whether quadrat-level trait distributions shift in a predictable direction along the soil fertility and acidity gradients within each forest plot. This analysis was then repeated using torus translation simulations [29]. The procedure included two steps: 1) we moved the true soil map by 20-m increments two-dimensionally, but kept the above trees map still; 2) We recalculated the Pearson correlation between T Q and F Q based on 20620 m quadrats for each simulation and compared the observed and simulated correlation coefficients. If the rank of the r-value from the true quadrat was higher than 97.5% or lower than 2.5% of the ranks of the simulated r-values (two-tailed test), it was considered that T Q and F Q was significantly correlated. The torus translations maintained the observed spatial distribution of soil fertility and tree distributions, but break their observed dependence by shifting the observed soil fertility distribution on a torus relative to fixed tree distributions.
We performed the three correlation analyses for five plant traits (SM, LA, SLA, WD, H max ) and the soil properties of the first two principal component axes (see Results: Soil Properties). We therefore used the false discovery rate (FDR) approach to adjust p-values [30,31]. Except the torus translation simulation, all other reported p-values refer to the adjusted p-values. All analyses were performed in R 2.13.0 (R core team, 2011).

Soil Properties
The significant soil PCA axes combined to explain more than 70% the variation of soil nutrients at each site (Table 1). For GTS, the first axis explained 42.2% of the overall variation and represented a general soil fertility axis with negative loadings on most key limiting elements. The second axis explained 17.8% of the overall variation and represented an acidity index with a relatively large negative loading on pH (20.461) and positive loadings on Fe (0.579) and B (0.610). The third axis explained 12% of the overall variation (Table 1).
For BCI, the PC1 axis explained 55.3% of the overall variation. Similar to the GTS analysis, this axis was also generally indicative of a soil fertility index with most elements decreasing. Again similar to GTS, the PC2 axis for BCI explained more than 12% of the overall variation and was an acidity index with low pH (20.470) and N (20.595) and high Al (0.348) and Fe (0.336). The PC3 axis explained about 11% of the overall variation and captured a large correlation between Al (0.600) and P (0.692) ( Table 1). This is also similar to the third axis at GTS, which had large loadings of the same sign for Al and P (Table 1).
For both plots, the PC1 axis was negatively related to soil fertility as shown by the negative loadings of most essential nutrients. Therefore, a negative correlation between a mean trait value and the soil fertility PC1 axis meant that trait values were larger on more fertile soils. The PC2 axis was positively related to soil acidity as shown by the negative loading of pH on PC2 (or lower pH for larger values of PC2) (Table 1). Therefore, a negative correlation with the soil acidity PC2 axis meant that larger trait values occurred on less acidic soils. Given the similarities between the loadings of nutrients on the first two PCA axes in both plots and because of their interpretability as general fertility and acidity axes, the following will focus on these first two axes.

Species-level Relationships between Trait Values and Soil Properties
Species functional trait values (T S ) were unrelated to calculated mean species soil properties (F S ) at BCI ( Table 2). For GTS, after the false discovery rate (FDR) correction to p-values, LA was negatively related to the soil fertility axis PC1 (r = 20.220; p = 0.008) and LA (r = 20.275; p,0.001), SLA (r = 20.221; p = 0.008) and WD (r = 0.269; p,0.001) were significantly related to the soil acidity axis PC2 (Table 2). We repeated these analyses for individual soil variables (see Table S1& Table S2). In sum, in both plots several significant relationships were uncovered between individual soil variables and all traits except for H max at GTS and SM at BCI.

Phylogenetically Independent Species-level Relationships between Trait Values and Soil Properties
A second set of species level analyses were performed using phylogenetically independent contrasts (PICs) to account for the evolutionary non-independence of species. In the GTS plot, five significant relationships were found between trait contrasts and soil contrasts (F S from eqn. 1) ( Table 3). In particular, significant positive relationships were found between SLA and the soil fertility axis PC1 (r = 0.302; p,0.001), SM and the soil acidity axis PC2 (r = 0.197; p = 0.028) and WD and PC2 (r = 0.286; p,0.001). Negative relationships were found between LA and the soil fertility axis PC1 (r = 20.176; p = 0.036) and SLA and the soil acidity axis PC2 (r = 20.245; p = 0.007).
In the BCI plot, five significant relationships were also found between trait contrasts and soil contrasts ( Table 3). The soil fertility axis PC1 was positively correlated with SM (r = 0.352; p,0.001) and WD (r = 0.167; p = 0.010). The soil acidity axis PC2 was positively correlated with LA (r = 0.268; p,0.001) and SM (r = 0.274; p,0.001) and negatively correlated with WD (r = 20.202; p = 0.003). Similar to the non-phylogenetic analyses, we repeated all analyses using individual soil nutrients (see Table S3 & Table  S4). Similar to the above non-phylogenetic results, several significant relationships were uncovered in both plots between individual soil variables and all traits except for SM at GTS.

Relationships between Quadrat-Level Mean Trait Values and Soil Properties
A second goal of this study was to test whether our predictions regarding species-level trait -soil relationships scale-up to the quadrat-level. Fifteen of the 20 possible relationships were significant after the FDR correction to probability levels ( Table 4). In the GTS plot, there was a negative relationship between the soil fertility axis PC1 and LA ( values could be found in Figure S3. In the BCI plot, the soil fertility axis PC1 was negatively correlated with H max (r = 20.283; p,0.001) and positively correlated with SLA (r = 0.258; p,0.001) and SM (r = 0.113; p,0.001).The soil acidity axis PC2 axis was positively correlated with LA (r = 0.242; p,0.001) and H max (r = 0.057, p = 0.037), but negatively correlated with WD (r = 20.206; p,0.001). The spatial distribution of all traits and soil PC1 and PC2 values are shown in Figure S4. As with the species-level analyses, all analyses were conducted on individual soil nutrients (see Table S5 & Table S6). At the quadrat-level all traits were correlated with at least one individual soil variable in each forest plot.

Torus Translation Simulations
As there is substantial spatial auto-correlation in species distributions and soil nutrient levels, we re-analyzed all of the quadrat-level trait-soil relationships using a torus translation approach. In Table 5, we provide the rank of observed Pearson r-values in the distribution of the randomized r-values for the 10 predictions for both BCI and GTS. The rank value could be used to calculate the significance of the observed correlations. In particular, low ranks or p-values indicated stronger than expected negative correlation and high ranks or p-values indicated a stronger than expected positive correlation. In the GTS plot, the observed significant Pearson correlation r-values between LA and SM and the soil fertility axis PC1 were still significant in the torus simulation (see Table 5; p.0.975 and p = 1.000). The observed relationship between the soil acidity axis PC2 and LA and WD were also significant (see Table 5; p = 1.000 and p,0.025). In the BCI plot, none of the observed r-values were significant after implementing the torus translation simulations and the false discovery rate (FDR) correction to p-values (Table 5).

Discussion
The distribution of plant species and communities along broadscale environmental gradients is expected to be determined by the sorting of species along these gradients on the basis of their function (e.g. [4,[8][9][10][11][12]32,33]). If a similar sorting of species by function occurs on local scales, then this may explain the substantial level of interspecific variation within local sites [13]. An expected mechanism underlying these trends is that species with 'fast' leaf, seed and wood economies that have faster resource acquisition and demographic rates should prefer resource rich ends of the gradient and species with 'slow' economies that have slower resource acquisition and demographic rates should prefer the resource poor ends of the gradient. The present analysis tested these mechanistic predictions in two forest plots. While the distribution of some plant traits showed significant relationships with local soil gradients, the majority of the expected relationships were not supported. In the following we discuss the results in detail for each forest plot.

Relationships between Functional Traits and Soil Resource Axes in the Subtropical Gutianshan (GTS) forest plot
Our species-level correlation analyses of the Gutianshan (GTS) forest plot in China supported three of our ten predictions when considering the results of both the species-level analyses and the phylogenetically independent contrasts (PICs) ( Table 6). Specifically, leaf area was positively correlated with soil fertility, specific leaf area was negatively correlated with soil acidity (pH) and wood density was positively correlated with soil acidity (pH). Leaf area and specific leaf area (SLA) are known to be correlated with high rates of resource acquisition and growth [11,17,34,35]. For example, species with high SLA values have low structural investment and relatively high photosynthetic and respiration rates, whereas species with low SLA values tend to invest more on leaf structures and have relatively low photosynthetic and respiration rates [2,18,36]. It was therefore expected that plants  Figure S3 for the complete maps of other traits and the soil PC2 values for GTS plot and Figure S4 for maps of all traits and the soil PC1 and PC2 values for BCI plot. doi:10.1371/journal.pone.0034767.g002 Table 4. Pearson correlation coefficients between five functional traits and the calculated scores of the two significant principal components of soil fertility and acidity for both GTS and BCI at the quadrat-level (Leaf area and Seed mass of BCI plot are log 10 transformed).  with larger SLA were found in high nutrient soils, whereas the reverse was expected to occur in low nutrient-supply soil. The trade-off between wood density and species growth and mortality rates has been shown in previous studies [37][38][39][40]. Speices in shaded or arid sites gernerally have smaller vessels and thicker fiber walls, thereby increasing their wood density. In the GTS plot, light wooded species tended to be found on fertile soils suggesting that high resource environments favoured species that allocate less biomass per unit volume and that have higher growth and mortality rates.

Relationships between Functional Traits and Soil
Resource Axes in the Tropical Barro Colorado Island (BCI) forest plot The species-level correlation analyses of the Barro Colorado Island (BCI) forest plot in Panama found no support for our predictions regarding species traits and soil fertility or acidity in species-level analyses ( Table 6), while the PIC analyses provided support for all but three of our predictions. Thus we could only support the predictions that seed mass would be positively related to soil acidity and negatively related to soil fertility and wood density would be negatively related to soil fertility (Table 6). In the low resource environments, large seeds could provide more reserves for individuals early in their life cycle. Small seeds, on the other hand, have the potential advantage of greater dispersal ability and rapid growth in high resource environments [41][42][43][44].

Quadrat-Level Trait-Soil Relationships in the Two Forest Plots
A secondary goal of the present study was to determine whether our predictions regarding species-level trait relationships with soil fertility and acidity gradients would scale-up to the quadrat-level. Although we found many significant relationships, the majority of these were non-significant once we accounted for spatial autocorrelation (Table 6). For example, in the GTS forest plot, the relationships between most traits and soil fertility were significant, but after controlling for spatial autocorrelation via a torus translation analysis, only four were still significant and only three of our ten predictions were still supported. This finding was consistent with our findings at species-level. From this, we can infer that for these few trait-soil relationships, it may be that the observed local relationships scale-up to generate the regional scale relationships reported elsewhere.
The quadrat-level results from BCI yielded no support for our predictions once we accounted for spatial autocorrelation. The non-significant co-variation between traits and soil at the quadratlevel may be based on very weak relationships at species-level in the BCI forest plot. A possible explanation for the BCI results is that the location for this forest plot was chosen to be as homogeneous as possible and a large proportion of trees there occur in shaded environments [45]. Therefore, the most important factors influencing the sorting of plant traits may be light levels or other factors rather than soil nutrients. Thus, this level of homogeneity also highlights one weakness of our study. Specifically, while the forest plots being analyzed used standardize tree inventory protocols, they were not set up to standardize the level of environmental heterogeneity. Future comparative research into the relationship between traits and soil nutrient gradients should therefore seek to standardize the level of soil nutrient heterogeneity at the plot-level.
Ultimately, the relationship between traits and soil fertility might be moderated by additional environmental parameters not presently analyzed and likely by the difference in the breadth of various resource axes. This may explain the lack of strong traitenvironment relationships in this study. Besides, a potential reason for the different results for the two plots is the difference in terms of seasonality and associated harshness of the abiotic environment.
In summary, plant functional trait research has shown that plant traits vary predictably along broad-scale climatic and soil gradients. The present research predicted that this variation might be explained partly by local-scale soil fertility and acidity gradients. Although we found leaf area and wood density had a consistent and predictable relationship with soil fertility both at species and quadrat-level for GTS, we failed to find support for most predicted relationships between plant traits and soil fertility and acidity axes. The correlations were calculated between the soil PC axes and the species-level or quadrat-level trait value. The table depicts whether the prediction was significant and supported the prediction (+), was significant and did not support the prediction (2) or was non-significant (NS). Phylogenetically independent contrasts (PICs) and torus translation simulations were utilized to correct for evolutionary non-independence in the species-level analyses and spatial auto-correlation in the quadrat-level analyses respectively. LA: leaf area; SLA: specific leaf area; SM: seed mass; WD: wood density; H max : maximum height. doi:10.1371/journal.pone.0034767.t006 In particular, the limited evidence for species-level associations between traits and soil fertility and acidity failed to scale up to the quadrat-level for both the GTS and BCI plots. The general lack of support for the predictions at the BCI forest plot may be due to limited heterogeneity in soil nutrients in this particular forest, but the same cannot be said for the GTS plot as it is quite heterogeneous with rugged terrain. In both plots it is clear that soil nutrients are not the only determinant of plant trait distributions and alternative resource axes, such as light, will have to be considered in future work. Ultimately, while some of our predictions regarding local-scale trait distributions, soil fertility and soil acidity were supported other factors likely play a larger role in determining the large interspecific variation in trait values in the forests studied. Figure S1 The phylogenetic tree constructed using DNA barcodes of the species in the GTS plot (See details on tree construction in the text).

Supporting Information
(TIF) Figure S2 The phylogenetic tree constructed using DNA barcodes of the species in the BCI plot (See details on tree construction in the text). Table S1 Pearson correlation coefficients between five functional traits and 13 soil nutrients for the GTS plot at the species-level (leaf area and seed mass are log 10transformed). (DOCX)