Interspecific trait integration increases with environmental harshness: A case study along a metal toxicity gradient

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. © 2020 The Authors. Functional Ecology published by John Wiley & Sons Ltd on behalf of British Ecological Society 1Laboratoire d'Ecologie Végétale et Biogéochimie, Université Libre de Bruxelles, Bruxelles, Belgium 2Environmental Change Institute, School of Geography and the Environment, University of Oxford, Oxford, UK 3Department of Forest, Nature and Landscape, Biodiversity and Landscape Unit, University of Liège, Gembloux Agro-Bio Tech, Gembloux, Belgium 4Ecology, Restoration Ecology and Landscape Research Unit, Faculty of Agronomy, University of Lubumbashi, Lubumbashi, Democratic Republic of Congo


| INTRODUC TI ON
Plant functional traits (sensu, Violle et al., 2007) are increasingly used to unravel assembly rules in plant communities. The functional structure of a plant community is generally described in terms of its species' trait mean values and trait dispersion (e.g. Bernard-Verdier et al., 2012). Both are driven by biotic and abiotic filters often linked to dramatic changes in community structure along environmental gradients (Kraft & Ackerly, 2013). However, traits do not vary independently from each other. Interspecific trait integration (ITI, the multivariate covariation between traits among species) might reflect functional trade-offs and constrain the space of functional variation within a community (Tilman, 1988). In a way, it relates to the concept of phenotypic integration, used in evolutionary biology to describe the covariation between phenotypic traits within species (Pigliucci & Preston, 2004).
Covariation patterns among functional traits have been extensively explored across species at large spatial scales. For example, the 'Leaf Economics Spectrum' (LES) has become a keystone concept of plant functional ecology (Wright et al., 2004). All plant species can be ordinated along a continuum of sets of traits: from long-lived organs with slow resource acquisition, to short-lived, fast acquisitive ones. The 'Whole Plant Economic Spectrum' (Reich, 2014) suggests that trait values for leaves, stems and roots should be highly correlated with each other and with relative growth rate, being reducible to a single axis of trait variation. Recently, Díaz et al. (2016) found that two orthogonal functional dimensionsnamely size and resource use-could be used to characterize plant functional diversity. The interpretation of community assembly mechanisms often relies on traits of the 'Leaf-Height-Seed' (LHS) scheme (Westoby, 1998), based on the idea that these traits reflect different functional axes and are therefore weakly correlated (Laughlin, Leppert, Moore, & Sieg, 2010). The importance of the multidimensionality of the phenotype has already been emphasized (Laughlin, 2014) and is often considered when investigating community assembly and functioning (e.g. Blonder, 2018;Laliberté & Legendre, 2010).
Although trait covariations have been acknowledged to reflect fundamental functional trade-offs, constraining strategies at large spatial scales, it is still unclear to what extent these covariation patterns hold true along local environmental gradients (Messier, McGill, & Lechowicz, 2010;Wright & Sutton-Grier, 2012). At fine spatial scales, Messier, Lechowicz, McGill, Violle, and Enquist (2017) and Messier, Violle, Enquist, Lechowicz, and McGill (2018) showed that the classical relationships between traits were not conserved in tree saplings communities, either at the interspecific level or at the intrapopulation level. Instead, it was suggested that traits related to different functions can form a network of community-specific trait interactions. Furthermore, in tropical forests communities, functions such as leaf economics and hydraulics spectrums can be decoupled (Li et al., 2015), hence suggesting many viable trait combinations and a low ITI. Nevertheless, the degree of trait covariation observed at the global scale has also been hypothesized to vary along environmental gradients at fine spatial scales. Trait integration was reported to increase along environmental stress at the community (Chapin III, Autumn, & Pugnaire, 1993), the intraspecific levels (Gianoli, 2004;Waitt & Levin, 1993) or at both (Read, Moorhead, Swenson, Bailey, & Sanders, 2014). As the environmental constraint increases, and because of the limited amount of energy to be allocated to growth, survival and reproduction, plant phenotypes are increasingly constrained by functional trade-offs. Therefore, the multidimensionality of functional traits per se could also influence species coexistence (Kraft, Godoy, & Levine, 2015;Laughlin & Messier, 2015).
Along local environmental gradients, ITI has been little investigated in comparison to trait mean values and other indices of functional diversity (see Mouchet, Villéger, Mason, & Mouillot, 2010).
Recently, Dwyer and Laughlin (2017) measured the change in trait covariation along an aridity gradient to quantify the effect of environmental filters on functional trade-offs. They showed that the positive trait covariation between vegetative height and seed mass within communities was a better predictor of species richness than other univariate or multivariate indices of functional diversity.
Consequently, trait covariation patterns may carry important information about the functional strategies selected along environmental gradient both within and between species and the relevance of ITI to understand the functional assembly of communities should be further investigated in a variety of environments.
In this paper, we explore the variation in ITI in metallophyte (i.e. species growing on metal-rich soils) communities along a metal toxicity gradient in the south of D.R. Congo. These savannas are unique from the floristic and functional points of view: a broad gradient of soil Cu and Co concentration drives a species and life forms turn-over in the communities (Boisson et al., 2020;Delhaye et al., 2016). This results in a decrease in size-related trait values, a change to fast resource acquisition strategies, and high leaf metal contents at the community level, while functional diversity decrease along with increasing soil Cu content . Furthermore, some species have a very narrow niche and are found strictly on soil with high or low Cu content, while others present a broad niche and occur along most of the gradient (Boisson et al., 2020). Yet, very little is known about global trait covariation patterns in the vegetation of metalliferous environments which could provide valuable information on plants' strategies to cope with the excess of metal (Faucon et al., 2016;Lange et al., 2017). Here, we analyse the covariation of functional traits, in the entire plant species pool of a metalliferous ecosystem and in local communities along a Cu content gradient.
We use bivariate and multivariate indices of ITI to ask two specific questions: 1. How do the traits of metallophyte species compare in the general whole plant economics spectrum? We expect to find the well-known relationships observed between traits at the global scale, such as the LES. We also expect a trade-off between plant metal content and other functional axes such as resource acquisition or conservation.
2. Is there an increase in ITI along the gradient of increasing toxicity that could explain the decrease in species richness observed? We expect an increase in functional trade-offs with the increasing soil metal content. We also expect the strong negative and nonlinear relationship between trait integration and species richness predicted by theory (Dwyer & Laughlin, 2017).

| Sampling
The study site is located on the north part of Fungurume V hill (S 10°37ʹ00ʺ E 26°17ʹ20ʺ), a 3.3 ha site situated approximately 200 km WNW of Lubumbashi in the D.R. Congo (Séleck et al., 2013).
The climate of the region is subtropical with a marked dry season (May-September) and with annual rainfall of 1,300 mm concentrated during the rainy season (November-April). The hill displays a strong gradient with soil Cu and Co concentrations increasing from the base to the top (92-6,736 mg/kg Cu and 10-927 mg/kg Co, respectively; acetate-EDTA extractible concentrations). Both concentrations are highly correlated (r = 0.92), representing a single, strong gradient of metal toxicity . Community composition and soil data were extracted from Séleck et al. (2013). Briefly, eighty-four 1-m 2 plots were installed, separated from each other by an elevation of 5 m, along six parallel transects from the base to the top of the hill. In each plot, the relative soil cover of each species was evaluated, and a composite soil sample was collected. In the present study, we use the acetate-EDTA extractible concentration of soil Cu as the gradient along which the ITI is studied. Reasons for this are that soil Cu contributed the most to the first axis of a PCA of all soil variables, and it constitutes a steep gradient of metal toxicity and rockiness , which, in turn, drives the composition of communities (Séleck et al., 2013) and species turn-over (Boisson et al., 2020).
For the analyses, we ordered the 84 plots according to their soil Cu content and dissected the gradient into 21 groups of four plots each. The resulting 21 composite communities prevent monodominance of large-stature species or local extreme stochasticity in species composition due to the size of the plots. In all communities, all traits were available for at least 18 species (see further).
Soil Cu content value of each composite community was calculated as the mean of the four plots. For functional trait measurements, we selected 65 abundant and representative species of the copper flora from the Katangan Copperbelt (see http://www.coppe rflora. org). This sample encompasses all dominant species, accounting for more than 90% of the total cover in each plot (Pakeman & Quested, 2007) but also some rare species, as these have been shown to present original and complementary functional trait combinations . Eight functional traits (vegetative height, leaf area, rooting depth, seed mass, specific leaf area, leaf dry matter, leaf Cu and Co contents) were measured for each species on 5-30 individuals distributed along the whole species-specific range of the Cu gradient (Boisson et al., 2020), following standardized protocols (Pérez-Harguindeguy et al., 2013 and see Appendix S1 for detailed trait measurement procedure). These traits (Table 1) encompassed the LHS strategy scheme (Westoby, 1998) as well as the LES (Wright et al., 2004), and leaf metal content related to the metal tolerance strategy (Baker, 1981). Traits values were averaged at the species level as intraspecific variation in this ecosystem has been showed to be low compared to interspecific differences (Delhaye et al., 2016).

| Statistical analyses
In order to achieve a normal distribution, LA, SLA, SM, Cu leaf and Co leaf were log transformed while VH and LDMC were square root transformed. All traits were standardized to a mean of 0 and a unit standard deviation prior to analysis. We investigated the relationships between traits using a principal component analysis (PCA) of the averaged trait values per species. To investigate the bivariate relationships between traits within each community, we computed a bivariate ITI index. For each pair of traits and in each community, a Spearman rank correlation coefficient was calculated. To avoid any potential bias related to the different numbers of species among communities, each coefficient was calculated on 1,000 subsamples of 18 species (i.e. lowest species richness among the local communities), weighting the probability to select a species by its relative abundance in the local community, and was then averaged for each community.

TA B L E 1 Description of the functional traits averaged at the species level, with minimum (Min) and maximum (Max) values
We tested for monotonous trends of the bivariate ITI index along the soil Cu content using Spearman rank correlation. A Šidák correction was used to adjust p values for multiple tests (Šidák, 1967), changing the significance threshold from 0.05 to 0.002. All analyses were carried out using the R software (R Core Team, 2019).

| Interspecific traits integration
For the PCA including all traits and species, the first three principal components (PC) accounted for 38%, 20% and 12% of the dataset total variance (Figure 1; Table S1 for loadings on axes and eigenvalues). The first PC was highly positively correlated with Cu leaf and negatively correlated with SM and VH. The second PC was highly positively correlated with LDMC and the third PC was positively correlated to SLA and Co leaf .
Correlations between traits were generally weak (see Table S2).
However, a few traits were consistently correlated within the whole  Table 1 for abbreviations and Table S1 for loadings Cu leaf Co leaf rooting depth were positively correlated to one another (ρ = 0.41) and to seed mass (ρ = 0.54 and 0.46) while negatively correlated to SLA (ρ = −0.4). Cu leaf and Co leaf were also significantly correlated (ρ = 0.52), reflecting the simultaneous tolerance of high foliar levels for both metals. Interestingly, we found no significant correlations between traits expected to be related to a same function such as plant size (plant height and leaf area or rooting depth) or LES (SLA and LDMC: ρ = −0.04).

| Variation of traits integration along the Cu gradient
Both multivariate ITI indices highlighted a sharp increase in phenotypic integration along the Cu gradient (sesITI range ρ = 0.88, p < 0.001, Figure 2a; sesITI sd ρ = 0.78, p < 0.001, Figure S1). In com- There was an increasingly negative association between vegetative height and SLA and between leaf area and SLA and Cu leaf .
Two traits associations (seed mass with leaf area and LDMC) decreased as soil Cu content increased (correlation coefficient getting closer to zero as soil Cu content increases). The correlation coefficient between Cu leaf and rooting depth changes from positive to negative, but with small, non-significant values.
Many trait associations were constant along the whole gradient, confirming the correlations observed within the whole species pool.
Among the traits showing the strongest association, Cu leaf and Co leaf were highly positively correlated, and both negatively correlated with vegetative height. Rooting depth was positively correlated with leaf area and seed mass and negatively with SLA. Finally, some trait F I G U R E 2 (a, c, e) Variation of the three multivariate indices (sesITI range , sesFRic and sesFDis) along the Cu gradient. (b, d, f) Association between species richness and the same indices. Each point represents one local community. The sesITI range index is the standardized effect size of the multivariate ITI calculated as the range of the eigenvalues of the PCA on eight traits in each community. sesFRic and sesFDis calculations are detailed in Appendix S2. The traits used are vegetative height, leaf area, rooting depth, leaf dry matter content, specific leaf area, seed mass, leaf Cu and Co content; see Table 1  associations were close to zero along the whole gradient (e.g. LDMC and SLA or Co leaf and rooting depth, see Figure 3).

| D ISCUSS I ON
The variation in covariation patterns between functional traits within communities and along abiotic gradients was recently suggested to be an important mechanism explaining species' distribution and taxonomic richness in communities. Here, this hypothesis was addressed using a tropical savanna developed on a steep gradient of Cu toxicity. We showed that, for a set of highly specialized species, some traits related to different functions are correlated. Further, the environmental filtering effect of the increasing Cu toxicity resulted in a strong increase in the ITI, increasing the association between some groups of traits that are not especially related to the same function. The results also supported the theoretical expectation that an increased trait covariation would drive a decrease in species richness (Dwyer & Laughlin, 2017). It is also worth mentioning that, regarding the bivariate ITI, while some trait associations were conserved along the whole gradient; other traits appeared to be consistently unrelated.

F I G U R E 3
Variation of the degree of association (Spearman correlation index) between the trait pairs along the Cu gradient. Each dot is the correlation coefficient of the two corresponding traits among all species of one community. The variation of the degree of association between traits along the Cu gradient is tested with Spearman correlation coefficient. Significant linear trends are illustrated by a black plain line (adjusted p value for multiple tests: p < 0.002)

| Functional traits integration in metallophytes
Our results did not support the idea of a tight covariation of all trait dimensions (Reich, 2014) nor the two orthogonal axes related to size and resources use (Díaz et al., 2016) observed at global scale.
We did not find strong support for the LHS strategy scheme, as smaller plants had larger SLA values. Furthermore, seed mass was highly correlated to leaf area and rooting depth, the latter traits reflecting plant size and resource capture axes. Instead, at the scale of a local plant community, we showed that traits related to different functions, such as size, resource use and metal tolerance, were moderately correlated, as suggested recently by . Specific leaf area and LDMC were not consistently related along the whole gradient suggesting that they were not related to the same function in our community, contrary to expectations (Laughlin et al., 2010;Wilson, Thompson, & Hodgson, 1999). This could be explained by the fact that LES traits are considered to reflect differences in leaf life span (Funk & Cornwell, 2013;, while most species in this study presented a similar leaf life span, due to the strong constraints imposed by the marked seasonality. This might also reflect the fact that LDMC is more strongly related to a leaf hydraulics axis, which was suggested to be decoupled from a leaf photosynthesis axis, here represented by SLA (Li et al., 2015). We did not find a strong positive correlation between seed mass and plant height as previously suggested (Cornelissen, 1999;Hodgson et al., 2017).
Larger-seed species, however, did have deeper roots and larger leaf areas, suggesting cohesive variation in some resource acquisition abilities across species. In these annually burned savannas, where biomass removal reduces competition for light (Zimmermann, Higgins, Grimm, Hoffmann, & Linstädter, 2010), competition could be more related to leaf area and rooting depth, advocating for the measure of several traits expected to be related to a same function (e.g. Aguirre-Gutiérrez et al., 2019).
Leaf metal concentrations are related to the strategy of metal tolerance (Baker, 1981). They can also be involved in the defence against herbivory (Behmer et al., 2005;Pollard & Baker, 1997) and could therefore be subject to trade-offs with other functions (Ernst, Verkleij, & Schat, 1992). When considering the entire species pool, smaller plants had a larger leaf Cu content. This suggests a tradeoff between interspecific competitive ability and metal accumulation (Ernst et al., 1992); the sequestration of metal in the plant could incur a cost that would not allow larger species to reach reproductive maturity in a short growing season. However, this result could also support the idea that leaf metal accumulation is a mechanism of biotic interaction (allelopathy) that would benefit poorly competitive, small species (Mohiley, Tielbörger, Seifan, & Gruntman, 2020). There were only moderate interspecific correlations between leaf metal concentrations and LES traits, potentially indicating that the physiological mechanisms required for metal tolerance within the leaves do not interact with resources use. However, these relationships should be investigated more thoroughly in other communities developing on soils with high levels of metal toxicity, which would benefit from the development of databases of traits (e.g. TRY) and of plant metal concentrations (Reeves et al., 2018). This should also be clarified by investigating experimentally physiological costs and effect on key plant functions of the different metal tolerance strategies.

| Variation in trait integration along the gradient
The central result of our study is that the strength of the relationships between some traits depends upon the abiotic environment and results in higher functional integration in harsher environments.
This supports the idea that environmental harshness reduces the quantity of viable trait combinations in the community. The environmental filter would therefore not only act on mean trait values but also on ITI. In turn, as suggested by Dwyer and Laughlin (2017), increasingly stronger filtering effects on ITI could account for the gradual decrease in species richness observed along the soil Cu content gradient Séleck et al., 2013). This filter effect would explain the decreasing proportion of species displaying a viable combination of traits with increasing Cu toxicity. Physical factors have been previously recognized to drive functional tradeoffs within communities by imposing energetic constraints between physiological functions, resulting in a more integrated phenotype (Boucher, Thuiller, Arnoldi, Albert, & Lavergne, 2013;Webb, Hoeting, Ames, Pyne, & LeRoy Poff, 2010). This also possibly provides a mechanism explaining the lack of trait coordination observed within certain communities (e.g. Silva, Souza, Caliman, Voigt, & Lichston, 2018).
This variation in multivariate ITI is the result of an increase in the bivariate trait correlations along the gradient, although this is not true for all trait pairs. This further supports the hypothesis that multivariate analyses can hide subtle patterns in traitenvironment associations (Butterfield & Suding, 2013). The bivariate trait association analyses also bring some complementary insights into the trait covariation patterns. For example, the positive correlation between SLA and leaf Cu content observed among the species is not constant along the whole gradient: the association between the two traits increases with soil metal concentration.
This reflects an increase in the proportion of metal rich, fast growing annual species on Cu-rich soils . We could hypothesize that accumulating and storing metals in the leaves is the most efficient strategy to ensure a rapid growth of annuals, while perennial species-with slower growth and flowering only after several years-might invest in long term, exclusion strategies.
In this sense, it was demonstrated that metallicolous population of Noccaea caerulescens, a hyperaccumulator species, present a shorter life cycle than non-metallicolous populations containing low leaf metal content (Dechamps et al., 2011). Furthermore, the associations between traits related to different functions in the LHS scheme (e.g. SLA and size) are also dependent on the environmental harshness, suggesting the necessity to reassess large patterns of traits covariation (e.g. Díaz et al., 2016) in relation to local environmental conditions. Interestingly, contrary to Dwyer and Laughlin (2017), we did not find any increase in the association between seed mass and vegetative height with the environmental harshness, but only a weak, non-significant, positive correlation along the whole gradient. This expected association was hypothesized to be due to trade-off between fecundity, germination and survival (Dwyer & Laughlin, 2017).
This result is likely due to the potentially weak contribution of competition for light in the interspecific interactions in these savannas.
The absence of this biotic filter in the communities may have relaxed the constraints on the association between these two traits related to competition for light. Further analyses could include the relationship between the reproductive output and the vegetative biomass as this ratio could be more informative than vegetative parameters to assess plant tolerance to metals (Dechamps, Lefèbvre, Noret, & Meerts, 2007).

| Comparing ITI and functional diversity indices
We found that both FRic and FDis values decrease with increasing soil metal concentrations, which is commonly interpreted as the effect of an environmental filter (Mason et al., 2013). However, these indices seem to perform poorly in predicting species richness variation among communities, as reported previously for FDis (Dwyer & Laughlin, 2017). This poor predictive capacity appears to be particularly strong for FRic (Figure 2). The new approach presented here introduces the strength of multivariate functional trade-offs as a driver of species richness in communities. This approach is radically different from commonly used functional diversity indices such as those used in our comparison. Indeed, these functional diversity indices focus on the shape of the multivariate trait space: FRic is the hypervolume of trait space, as defined by Cornwell, Schwilk, and Ackerly (2006), while FDis is the mean distance of each species to the centroid of the trait space, taking species abundances into account (Laliberté & Legendre, 2010;Villéger et al., 2008). In this regard, both indices are multivariate and contain a component of trait covariance. However, we expect ITI to be independent from these indices, as several patterns of trait covariation can produce similar hypervolumes and distances from the centroid. Multivariate ITI explains the decrease in species richness observed along our toxicity gradient better than functional diversity indices. This suggests that environmental filters do not select trait values independently but rather selects particular trait combinations, expressed at the community level by a stronger ITI between some traits. Analysing ITI therefore represents a complementary approach to investigate the community trait space and the ecological processes affecting it. We suggest that rigorous tests of the redundancy/complementarity of functional trait-related indices should be conducted in order to integrate functional trait integration in the larger framework aimed at characterizing the functional diversity of communities and predict community responses to rapidly changing environments towards physiologically stressing conditions (e.g. climate change).

| CON CLUS IONS
The analysis of ITI is a promising tool to investigate variation in multidimensional functional niche of species along gradients of environmental constraints. Since selection acts on multidimensional phenotypes, species traits are not independent from one another. Although trait values and ranges within communities can be expected to vary along gradients of environmental conditions, trait integration is likely to be a crucial element to uncover the ecological processes driving community assembly. We found that well-known relationships between traits at global scale may not be conserved in local communities, and more so in heterogeneous environments. These findings suggest that trait-trait relationships should be investigated in relation with environmental variables both at the community level and along larger environmental gradients. This could further provide some mechanistic explanations for the discrepancies observed between some interspecific and intraspecific trait covariation patterns (e.g. Rosas et al., 2019).