Diversity in form and function: Vertical distribution of soil fauna mediates multidimensional trait variation

It has been widely recognized that species show extensive variation in form and function. Based on species' attributes, they can be positioned along major axes of variation, which are often defined by life-history traits, such as number of offspring, age at maturity or generation time. Less emphasis has been given in this respect to tolerance traits, especially to resistance to abiotic stress conditions, which often determine community (dis)assembly and distribution. Soil fauna species distribution is governed to a large extent by environmental conditions that filter communities according to functional traits, such as abiotic stress tolerance, morphology and body size. Trait-based approaches have been successfully used to predict soil biota responses to abiotic stress. It remains unclear, though, how these traits relate to life-history traits that determine individual performance, that is, reproduction and survival. Here, we analyse patterns in multidimensional trait distribution of dominant groups of soil fauna, that is, Isopoda, Gastropoda and Collembola, known to be important to the functioning of ecosystems. We compiled trait information from existing literature, trait databases and supplementary measurements. We looked for common patterns in major axes of trait variation and tested if vertical distribution of species in the soil explained trait variation based on three components of trait diversity (trait richness, evenness and divergence). Our results showed that two to three axes of variation structured the trait space of life-history and tolerance traits in each of the taxonomic groups and that vertical distribution in soil explained the main axis of trait variation. We also found evidence of environmental filtering on soil fauna along the vertical soil distribution, with lower trait richness and trait divergence in soil-dwelling than in surface-living species. Our study was partially limited by the lack of detailed trait measurements for the selected taxonomic groups. In this regard, there is an urgent need for standardized trait databases across invertebrate groups to improve trait-based diversity analysis and fill gaps in the mechanistic understanding behind trait distribution, trait filtering and the link with species fitness and performance.


| INTRODUC TI ON
Evolution has led to an astonishing biological diversity in the Earth's terrestrial ecosystems. A major component of biodiversity is the variation in morphological, physiological or phenological features of organisms, also defined by Violle et al. (2007) as functional traits, which impact fitness indirectly via their effect on growth, reproduction and survival. Plants, animals as well as micro-organisms vary greatly in their allocation of resources to growth, survival and reproduction, giving rise to a wide variety of form and function, even on a small spatial scale. Understanding the patterns of diversity and variation, based on determinants such as phylogenetic history, geographic position, dispersal ability and habitat characteristics, is important for our ability to make predictions of species distributions, community composition and ecosystem functioning under environmental change (Lavorel & Garnier, 2002;McGill, Enquist, Weiher, & Westoby, 2006).
Ecologists and evolutionary biologists alike have recognized for a long time that traits do not vary freely within and among species but, on the contrary, the variation of functional attributes within plants and animals occurs in integrated trait syndromes. In plants, the leaf economics spectrum is a good example of such universal syndrome of key chemical, structural and physiological properties describing a set of trade-offs among traits related to plant carbon, nitrogen and phosphorus balance (Wright et al., 2004), and resulting in predictable relationships between traits. In animals, reproductive, developmental, dispersal, synchronization and life-history traits form an integrated response to particular ecological problems Siepel, 1994;Verberk, Siepel, & Esselink, 2008).
For instance, short development time and high dispersal rate are coupled in freshwater macroinvertebrates in ephemeral habitats, whereas slow growth and high adult longevity are found in environments with constantly harsh conditions (Verberk et al., 2008). Trait syndromes are often characterized by trade-offs, which may result from the allocation of limited resources between key life-history traits, such as size and number of offspring, age at first reproduction and growth rate (Le Lann et al., 2014;Liefting, Grunsven, Morrissey, Timmermans, & Ellers, 2015;Reich et al., 2003;Stearns, 1989).
Covariance between traits may also result from pleiotropic effects, that is, when one gene influences two or more seemingly unrelated phenotypic traits (Stearns, 1989). For example, physiological adaptations to abiotic stress such as drought have antagonistic effects on tolerance to inundation (Dias et al., 2013). However, less is known on how tolerance traits relate to life-history traits that determine individual growth, reproduction and survival.
Defining the major dimensions of trait space is a fruitful way to get a better understanding as to how ecological conditions shape the evolutionary trajectories of species and their range distributions ). An early paradigm describing such division of lifehistory strategies is the r/K selection theory (MacArthur & Wilson, 1967;Pianka, 1970). A similar axis of variation along which traits vary is the fast-slow continuum of reproductive traits, with species that mature early, have large reproductive rates and short generation times occupying the "fast" end of the continuum and those with the opposite suite of traits occupying the "slow" end (Blackburn, 1991;Franco & Silvertown, 1996;Promislow & Harvey, 1990;Read & Harvey, 1989;Southwood, 1988). In other studies, it was proposed that life-history variation can be characterized by two independent axes that typically inform on "the speed of life," one capturing variation in longevity and mortality schedules, and the other reflecting reproductive strategies in terms of timing of reproductive bouts (Bielby et al., 2007;for mammals;Salguero-Gómez et al., 2016;for plants). A successful characterization of the major axes is particularly helpful if it enables the use of few, easily measurable traits to represent species position along these axes. Plant ecologists have made important progress in explaining effects of climate change on species distribution, range shifts under environmental change and community (dis)assembly across ecological scales by applying this approach (Carreño-Rocabado et al., 2016;Cornwell & Ackerly, 2009;Garnier et al., 2004;Reich, 2014). However, particularly for major animal groups such as invertebrates, efforts to describe the dimensions of variation in form and function at community scale have only started recently (Fountain-Jones, Baker, & Jordan, 2015;Moretti et al., 2017).
Here, we aimed at characterizing major axes of variation in traits, including morphological, physiological, behavioural and life-history traits (sensu Moretti et al., 2017) for soil fauna. Soil fauna, such as earthworms, millipedes, isopods and springtails, are key to the functioning of soils (e.g., Bardgett & van der Putten, 2014). Hence, predicting the population performance of soil fauna in a changing environment is crucial to our understanding of ecosystem functioning and service provision. First, we compiled trait values from existing literature and trait databases or performed supplementary measurements for three groups of soil fauna: Isopoda, Gastropoda and Collembola. We then tested if the variation in trait space can be captured by a few main axes of variation using principal component analysis (PCA) and if the main axes correspond to environmental conditions that govern species distribution by filtering them according to abiotic tolerance. In soil, ecological conditions vary analysis and fill gaps in the mechanistic understanding behind trait distribution, trait filtering and the link with species fitness and performance.

K E Y W O R D S
Collembola, evolutionary trade-off, functional trait, Gastropoda, Isopoda, life history, soil trait diversity, vertical distribution steeply along a small-scale vertical stratification gradient (Berg & Bengtsson, 2007;Krab, Oorsprong, Berg, & Cornelissen, 2010), so that trait variation is expected along this gradient. Vertical distribution in soil captures the most relevant environmental stress factors for soil organisms, which are humidity and temperature. Therefore, we tested if the major axes of variation differentiate species according to their vertical distribution in soil, that is, the soil layer at which species live, and if this explains variation in tolerance traits of soil invertebrates to abiotic conditions. For Isopoda, a phylogeny of the most common northwest European species is available (Dias et al., 2013), which we used to assess phylogenetic signal in each trait and to perform a phylogenetically informed PCA taking into account non-independence of lineages; for the other groups, phylogenetic information was insufficient to do so.
In addition, to gain a more comprehensive understanding of changes in the multidimensional trait distribution of species along the vertical soil stratification gradient, we analysed three independent and complementary components of trait diversity: trait richness, evenness and divergence (Mason, Mouillot, Lee, & Wilson, 2005). Trait richness informs about the amount of functional space occupied by the species living in the three soil layers (e.g., Fontana, Petchey, & Pomati, 2016;Mason et al., 2005). Trait evenness quantifies how regularly distributed species are in the functional space defined by multiple traits (Fontana et al., 2016). Trait divergence measures the degree of trait dispersion around the centroid of the distribution (Laliberté & Legendre, 2010). These measures allowed us to investigate the diversity of ecological strategies potentially found in each soil layer across Europe, although we did not sample real communities in the natural environment. We hypothesized that specific environmental conditions can select for a reduced number of trait combinations (e.g., Cornwell, Schwilk, & Ackerly, 2006;Mouillot, Graham, Villéger, Mason, & Bellwood, 2013): we expected this environmental filtering (sensu Götzenberger et al., 2012) in deeper soil layers to result in reduced trait space coverage (i.e., lower trait richness) and convergence of species towards specific trait combinations (i.e., lower trait divergence and trait evenness).

| Trait data
We obtained traits for three groups of invertebrates that are commonly found in soils and that predominantly feed on detritus, that is, Isopoda, shelled Gastropoda and Collembola (Table 1). For each taxonomic group, we compiled a database with trait values (data available from the Dryad Digital Repository: https://doi.org/10.5061/ dryad.m6dn0g8). Most traits were obtained from existing databases, supplemented with data from the literature or with new measurements. We had to balance decisions on which traits to include for each taxonomic group, because the PCA (see below) required a completely filled trait matrix or only few missing values (Dray & Josse, 2015). Including more traits would increase the information on the shape of the axes, but trait values were often only available TA B L E 1 Description of functional traits of Isopoda, shelled Gastropoda and Collembola, as used in the analyses only data on total fecundity were available. Therefore, we inferred clutch size of these two species from their total fecundity using the equation obtained by correlating the reported values for total fecundity and clutch size for the other species in our Isopod dataset. This relationship is linear (clutch size = 0.5517*total fecundity + 4.7545) and very strong (R² = .89). Additionally, we measured clutch size for a number of species present in our reference collection, but without available literature data: Armadillidium pictum (n = 2), A. pulchellum (n = 4), Haplophthalmus mengei (n = 6), H. riparius (n = 5), Miktoniscus patience (n = 1) and Trichoniscoides sarsi (n = 1). The sample size for clutch size was low for a number of species (Table S4). To ascertain that this did not bias our results, we also analysed the major dimensions of trait space (see below) with only those species that had a sample size larger than n = 5. The results are shown in Table S3 and were nearly identical to the analysis with the full dataset. All trait values are deposited in the Dryad Digital Repository: https://doi. org/10.5061/dryad.m6dn0g8 , and the samples sizes are listed in the Supplementary Material (Table S4). we selected maximal shell size, survival of dry period (equivalent to drought resistance) and inundation tolerance (equivalent to inundation resistance), age at maturity, longevity and clutch size (number of eggs per clutch) ( Table 1). The database was also used to collect information on vertical stratification of species. The authors of the database acknowledge that knowledge about some of the relationships described in the database may be incomplete or imprecise in nature (Falkner et al., 2001). Therefore, all traits in the database were ordinal, and Falkner and co-authors describe the affinity of each species for the different categories with a fuzzy coding system: 0 (=no association), 1 (=minor association), 2 (=moderate association) or 3 (=maximum association). For our analysis, we transformed this clas- weeks = 3; months = 4) and inundation resistance (low = 1; moderate = 2; high = 3). Finally, vertical stratification in soil was derived from microsite data coded in the original database and expressed as a categorical trait with levels "soil-dwelling" (if "epigeon" < 3, or "epigeon" = 3, and "hypogeon" > 1), "mixed depth" (if "epigeon" = 3 and "hypogeon" = 1, or "epigeon" = 3 and "hypogeon" = 0, and "among/ under surface debris" = 3) and "surface-living" (if "epigeon" = 3 and "hypogeon" = 0, and "among/under surface debris" <3  (1994), Fjellberg (1998Fjellberg ( , 2007 and Bellinger, Christiansen, and Janssens (2017). We selected the following traits: maximum body size, temperature preference, moisture preference, thermal niche breadth and mode of reproduction, and we compiled data on vertical stratification in soil ( Table 1). Most of the selected traits, with the exception of body size and mode of reproduction, were proxies for the functional traits that determine species performance. These functional traits were only available for a very limited set of Collembola species, and therefore, it was not feasible to include them in the analysis. Temperature preference was used as a proxy for heat resistance and estimated from the global distribution maps of species (Bellinger et al., 2017). We assumed that boreal species are cold-adapted and

| Major dimensions of trait space
A PCA was performed for each taxonomic group. Tables S1a-c report the traits included for each of the taxonomic groups. The scores of the PC axes were extracted to test whether each axis differentiated between species with different vertical distribution in soil, using a two-sample t test between PC values of soil-dwelling species and surface-dwelling species (for Isopoda) or by performing a oneway ANOVA (for Gastropoda and Collembola).

| Trait diversity
Using the same traits as included in the PCA, we comprehensively in Collembola traits caused coplanarity issues during the calculation of TOP; therefore, convex hulls were estimated using joggled input (for all groups, for consistency), which allowed to handle this kind of data by slightly moving each input coordinate to guarantee simplicial facets (e.g., three-dimensional triangles). As traits were expressed in different units, we standardized them (M = 0, SD = 1) prior to calculation to make sure they all had the same weight in determining trait diversity indices. The different soil layers included a variable number of species (Tables S1a-c

| Phylogenetic signal in traits and phylogenetically informed PCA
Phylogenetic relatedness may cause non-independence of lineages due to shared history (de Bello et al., 2015). Phylogenetic methods take non-independence of lineages into account and should be preferred if a phylogeny is available. For the taxonomic groups used in this study, a suitable phylogeny was only available for Isopoda (Dias et al., 2013), based on the 18S gene. For the other groups, phylogenetic information was insufficient to do so. We are aware that single-gene phylogenetic trees are subjected to topological variation (Castresana, 2007); however, the provided isopod tree was highly concurrent with taxonomic classification and morphological traits.
Therefore, we are confident that no bias was introduced using this phylogenetic tree. We used this phylogeny to assess phylogenetic signal in each of the isopod traits. We calculated Pagel's λ, a scaling parameter for the extent to which correlations in traits reflect their shared evolutionary history, for each trait used in our analysis.
Pagel's λ values range from 0 (no correlation) to 1 (the correlation expected under Brownian motion). We have included Figures S2-S6 with the Isopod traits mapped on the phylogenetic tree using the R package ape and phytools. We then carried out a phylogenetically informed PCA (R Package Phytools, Revell (2009)) on the same traits as used in the traditional PCA to obtain the major axes of trait variation while accounting for non-independence of species.
All analyses were performed in R, version 3.3.3.

| Isopoda
PCA on body size, drought and inundation resistance, walking speed and clutch size showed that the first two PC axes captured the spectrum of traits of Isopoda adequately. PC1 and PC2 contained a cumulative proportion of 84.5% of the variation (Figure 1). PC1 explained 65.6% of variation and was positively correlated with body size, clutch size and drought resistance, and negatively with inundation resistance. All four of these traits loaded to a similar extent on PC1. PC2 explained 18.9% of the variation in trait distribution, and walking speed showed a strong positive association with this axis.
We then tested if a species' score on the first two PC axes could be predicted from the vertical distribution of the species in soil. PC1 revealed significant divergence between soil-dwelling and surfaceliving species (t = −11.38, df = 13.71, p < .001), with negative PC1 scores associated with soil-dwelling species, showing a decreased body size, smaller clutch size, lower drought resistance and an enhanced inundation resistance. No difference between soil-dwelling and surface-living species was instead found for PC2 (t = 0.53, df = 17.2, p = .60).
The role of phylogenetic relationships in determining the spectrum of isopod trait combinations was weak for body size and clutch size, as indicated by a low Pagel's lambda for these traits (Table 2). In contrast, lambda was 0.72 and 0.63 for inundation resistance and drought resistance, respectively, suggesting a stronger role of phylogenetic ancestry for those traits. To ensure phylogenetic independence in the assessment of the axes of trait variation, we additionally performed a phylogenetically informed PCA on the same five traits. The result showed PC1 to be highly correlated with inundation resistance and to explain 91.4% of the trait variation. PC2 explained an additional 5.4% and was strongly associated with drought resistance (Table S2). Therefore, the two tolerance traits defined the dimensions of trait space even if phylogenetic dependence was accounted for.

| Gastropoda
The trait space of Gastropoda was captured by three main axes of variation. PCA on maximal shell size, inundation tolerance, survival of dry period, clutch size, age at maturity and longevity showed that the first three PC axes contained a cumulative proportion of 83.3% of the variation (Figure 2

| Trait diversity
In all taxonomic groups, the three trait diversity indices significantly varied among soil layers, but the pattern differed between groups ( Figure 4).
In Isopoda, trait richness (TOP) and trait divergence (FDis) were significantly lower in soil-dwelling than in surface-living species of Collembola across soil depth differed from the other two groups, which was mainly due to a strong effect of mode of reproduction (only used for Collembola). When this trait was excluded from the analysis, all trait diversity indices significantly decreased with depth in the soil ( Figure S1).
F I G U R E 2 Alignment of Gastropoda traits with the first two principal component analysis (PCA) axes. Six traits (i.e., maximal shell size, age at maturity, longevity, clutch size, inundation tolerance and survival of dry period) were measured for 169 terrestrial snail species included in PCA. The three clusters indicate the vertical distribution of species with soil depth: soil-dwelling species (in red), mixed depth (green) and surface-living species (in blue). Note that vertical distribution was not a variable included in the PCA. Ellipses depict 95% confidence intervals around the mean, assuming a normal distribution. Therefore, ellipses that do not overlap are likely significantly different (α = 0.05) F I G U R E 3 Alignment of Collembola traits with the first and second principal component analysis (PCA) axes. Five traits (i.e., body size, mode of reproduction, moisture preference, temperature preference and thermal niche breadth) were measured for 278 springtail species included in PCA. The three clusters indicate the vertical distribution of species with soil depth: soil-dwelling species (in red), sub-surface-living species (green) and surface-living species (in blue). Note that vertical distribution was not a variable included in the PCA. Ellipses depict 95% confidence intervals around the mean, assuming a normal distribution. Therefore, ellipses that do not overlap are likely significantly different (α = 0.05)

| D ISCUSS I ON
The diversity in form and function of plants and animals can be characterized by major axes of life-history variation and reproductive strategies (e.g., Díaz et al., 2016). Here, we aimed at defining the major dimensions of trait space of soil invertebrates, a group that is impor- For the Isopoda, a suitable phylogeny was available, which allowed us to look for phylogenetic signal in the traits. The life-history traits clutch size and body size showed low phylogenetic signal, indicating that these traits were hardly phylogenetically constrained. This confirms earlier findings in plant life-history strategies where phylogenetic ancestry also played a minor role (Salguero-Gómez et al., 2016). However, for the tolerance traits inundation resistance and drought resistance phylogenetic relationship were more important in explaining trait variation. Other stress tolerance traits, such as upper thermal limit, have also been shown to contain a strong phylogenetic signal in insects (Araújo et al., 2013;Kellermann et al., 2012). Ignoring phylogeny would result in an incorrect assumption of statistical independence for those traits and potentially lead to elevated type I error in the PCA analysis (Revell, 2009 (Bielby et al., 2007;Díaz et al., 2016;Salguero-Gómez et al., 2016).
In the three taxonomic groups we analysed, however, the major dimension of variation was not a fast-slow continuum but a vertical distribution in soil axis, although, in the Gastropoda, we found indication for a minor fast-slow axis of variation in PC3. Several reasons can be invoked to explain lack of evidence for a fast-slow continuum.
First, the trait space we explored was more limited than in previous studies (e.g., Bielby et al., 2007;Blackburn, 1991 In contrast, our analysis for all three taxonomic groups included tolerance and resistance traits to abiotic stresses, that is, inundation tolerance, drought resistance and thermal niche breadth, which to our knowledge have not been included in trait spectrum studies before. Tolerance traits are expected to evolve under selection of the abiotic conditions a species experiences in its habitat (Kellermann et al., 2012;Van Dooremalen, Berg, & Ellers, 2013). For soil fauna, vertical stratification in soil is one of the main determinants of abiotic conditions as there are steep gradients of temperature and moisture conditions across soil profiles (Berg & Bengtsson, 2007;Krab et al., 2010). The strong and consistent signature of vertical stratification on the first major axis of the PCA probably resulted from including tolerance traits in our analysis. Importantly, this was also followed by variation in life-history traits, such as clutch size, age at maturity and longevity. An important focus for future studies is to investigate how tolerance traits are correlated with life-history traits, also across other plant and animal species. Major habitat type has been found to be a weak, but significant, predictor of plant species distribution in functional trait space (Salguero-Gómez et al., 2016), although no direct measurements of tolerance traits were used in their study. Although tolerance traits directly underlie species performance, a drawback is that they are more difficult to compile as they require labour-intensive measurements under standardized conditions (Moretti et al., 2017).
Additional understanding of the responses of species to abiotic and biotic conditions was provided by our integrated analysis of the multidimensional trait space. We showed that trait diversity differed significantly between surface, sub-surface and soil environments. In Gastropoda and Isopoda, the results suggested environmental filtering, with lower trait richness and trait divergence in the lower strata.
We refer here to environmental filtering in a loose way, simply meaning a contraction of the trait space in a given environment, as we are aware that both abiotic and biotic condition (i.e., competition and predation) can determine community assembly (Cadotte & Tucker, 2017 (Chahartaghi, Scheu, & Ruess, 2006).
Interestingly, according to Petersen (2002), soil-dwelling Collembola species are small, produce few but large eggs and reproduce throughout the year compared to surface-living species, suggesting a coupling between tolerance traits and life-history traits.
As hypothesized, trait evenness tended to decrease with increasing depth (Gastropoda and Collembola when excluding mode of reproduction, see Figure 4 and Figure S1, respectively  (Clark et al., 2011). It is important to note, though, that traits have been traditionally investigated mainly in terrestrial plants, whose traits are relatively easy to measure following available standardized procedures (Pérez-Harguindeguy et al., 2013). In the case of arthropods, standard protocols for the measurement of terrestrial invertebrate functional traits are now available (Moretti et al., 2017), and they will hopefully encourage the effort of measuring enough individuals to adequately investigate intraspecific trait variation in the future.
Trait-based ecology has developed into one of the principal fields that address macroecological and evolutionary questions.
Early studies have successfully compiled large online trait databases to address macroecological questions such as the effects of climate change and stressors on species distribution and ecosystem processes (e.g., Berg et al., 2010;Diamond, Frame, Martin, & Buckley, 2011;Dias et al., 2013;Moretti & Legg, 2009

ACK N OWLED G EM ENTS
We thank the editors for the invitation to contribute to this issue and for their feedback. Three anonymous reviewers provided valuable comments on an early version of the manuscript. J.E. was sup-

DATA ACCE SS I B I LIT Y
Data from this study are available online at the Dryad Digital Repository: https://doi.org/10.5061/dryad.m6dn0g8 .