Functional and Taxonomic Diversity of Collembola as Complementary Tools to Assess Land Use Effects on Soils Biodiversity

Collembola have been proposed for several decades as a good model organisms to survey soil biodiversity; but most of the studies focused on taxonomic endpoints. The main objectives of this study are to compare the effects of the different land uses, including urban and industrial land uses, while using both collembolan functional and taxonomic biodiversity approaches. We collected data on 3,056 samples of Collembola communities across 758 sites in various land uses throughout France. The types of land use considered included all types of human activity from forestry to urban, industrial, traffic, mining and military areas, agricultural grassland, arable land, vineyards and urban vegetable gardens. In order to study functional and taxonomic biodiversity, we used community-weighted means, functional indices, species richness and density. When looking at collembolan functional diversity, urban and industrial soils appear clearly less diversified than when considering the taxonomic diversity. We suspect here a functional homogenization effect commonly reported in the literature for various organisms in urban ecosystems. Our study provides range of values for different taxonomic and functional indices of Collembola communities in a wide land use classification across France.


INTRODUCTION
Soil biodiversity, which can represent up to 25% of the terrestrial biodiversity (Decaëns et al., 2006), is usually considered as a key component of ecosystem functioning, providing many ecosystem services. In this context, land use represents one of the main global anthropogenic parameters modifying soil physico-chemical characteristics and filtering soil communities. Thus, Joimel et al. (2016) illustrated a gradient of intensity of soil anthropized, from forest, managed grassland, Frontiers in Ecology and Evolution | www.frontiersin.org 1 July 2021 | Volume 9 | Article 630919 cultivated agriculture, orchards and vineyards, urban vegetable gardens to urban, industrial, traffic, mining and military areas, thanks to fertility and contamination. Among soil biodiversity, Collembola have been used as bioindicators in a number of national soil monitoring programs across Europe (George et al., 2017) and provide a tool for assessing the effects of land use on soil biodiversity (Vandewalle et al., 2010). They participate in many soil functions including litter decomposition (Cortet et al., 2003), nitrogen and carbon cycles (Filser, 2002), and can also impact soil microstructure (Rusek, 1998). Moreover, Collembola are extremely common and dense in the upper soil.
It has been demonstrated by Joimel et al. (2017) that, contrary to the conclusions obtained from physico-chemical indicators, urban soils do not show the lowest soil Collembola biodiversity, especially compared to arable and vineyard cultivation areas, which appear to be the most harmful land use. However, this study was based on a taxonomic diversity indices and show only one facet of biodiversity. The functional diversity of soil invertebrates can be assessed by functional trait-based approaches, which are often perceived as a good tool for monitoring soil fauna (Reis et al., 2016). Functional traits defined as morphological, physiological, phenological or behavioral features measured at the individual level (Pey et al., 2014), are relevant imprints of their adaptation to the environment and/or can affect ecosystem properties (Loreau et al., 2001;Diaz et al., 2004). Even, if it is difficult to directly relate soil fauna traits to services, it is possible to construct a list of key traits related to three main ecological functions that are at the core of soil suitability for the provisioning of ecosystem services: (i) soil formation; for example, concerning the aggregation process, the size of ant workers' mandibles could determine the size of aggregates they are able to move (Dostál et al., 2005), (ii) nutrient cycling; body size appears to be a key functional trait as larger animals have generally the strongest effects on decomposition (Rusek, 1998). Moreover, Peguero et al. (2019) found that taxonomic, phylogenetic, and functional richness of springtail communities alone were able to explain up to 30% of the variation in annual decomposition rates, and (iii) pest control; Brousseau et al. (2018) demonstrated a better predictive power with functional traits than the commonly used predator/prey body size ratio.
Recently, several studies have been conducted on European and national scales in order to assess the effect of land use or soil managements on soil functional diversity (Salmon and Ponge, 2012;Salmon et al., 2014;Martins da Silva et al., 2015), although these efforts have mainly concerned non-urban soils and often less or moderately anthropized soils (forest, agriculture, and permanent crops). With only 1% of soil studies published dealing with urban soils (Guilland et al., 2018), a considerable knowledge deficit precludes obtaining reference values for soils (Guilland et al., 2018). Such a knowledge gap leads to an insufficient consideration of urban soils by city managers and urban planners (Blanchart et al., 2017). Soils from urbanized areas are defined as SUITMAs, because they include soils of urban, sensu stricto, but also industrial, traffic, mining, and military areas. This definition refers to a large number of soil types of strongly anthropized areas . One of the main questions related to SUITMAs concerns biotic homogenization (McKinney and Lockwood, 1999). A previous study concerning soil Collembola, conducted in urban vegetable gardens, demonstrated that, contrary to plants, soil taxonomic, and functional Collembola communities remained quite dissimilar between three cities in France (Joimel et al., 2019). However, these homogenization patterns observed must now be compared to other land uses to obtain a global picture.
A large Collembola community dataset has been collected over 20 years from 13 different research programs concerning various land uses (forest, agricultural, urban and industrial) in France. Our main objectives are to compare the effect of the different land use, on both Collembola functional and taxonomic biodiversity, including SUITMAs. We hypothesize that functional homogenization is greater in urban systems compared to other land uses, even when the species composition changes, as a consequence of the selection of traits associated with high adaptive value in the urban environment (convergence of traits).

Dataset Collection
A dataset of 3,096 samples from 758 sites was obtained from several research programs (Table 1). To avoid bias, data was obtained from various projects using the same extraction and sampling methods. Microarthropods were extracted over a period of a week from intact soil cores (5 cm depth, 6 cm diameter), using a high-gradient Macfadyen extractor (ISO, 2006). Between 3 and 6 soil cores (soil samples) were sampled from each site (replicate), according to each sampling strategy and the data collected was averaged. For the present study, the data pertain exclusively to collembolan communities in topsoils to a depth of 5 cm (more information in Joimel et al., 2017). A total of 143 taxa of Collembola were identified to the species level. Values of density (abundance of each taxon per m 2 ) and taxa richness were recalculated from Joimel et al. (2017).
The land use examined consists of the following types: forest, grassland, arable land, vineyard, urban vegetable garden, and SUITMA (Soils of urban, industrial, traffic, military, and mining areas). Likewise, as in the study of soil chemical properties (Joimel et al., 2016), the SUITMA group comprises soils used by the mining industry, solid, and liquid waste dumping, as well as habitation and road construction. In contrast, urban garden soils were excluded from the SUITMA category and constitute a separate type of land use in the present study, as they represent specific geochemical characteristics ( Table 2).

Functional Trait Composition
Functional traits were selected according to hypotheses linking land use characteristics and adaptation to subhabitats (soil or soil surface). We focused on traits relating to i) vertical distribution of Collembola within the soil profile (Salmon et al., 2014), i.e., body shape, presence of visual organs (ocelli and post-antennal organ), the presence of functional furcula or scales and pigmentation, Frontiers in Ecology and Evolution | www.frontiersin.org and ii) their potential resistance to disturbance or stress caused by urban environment, i.e., parthenogenetic reproduction.
Trait values were extracted from the BETSI database (Biological and Ecological functional Traits of Soil Invertebrates) 1 . All data are accessible on our institutional data repository: https://doi.org/10.15454/UU2FQT. Traits were split into attributes by a fuzzy coding approach (Pey et al., 2014). As recommended by Bonfanti et al. (2018), we selected only Pan-European keys in order to limit intraspecific variability. Body shape, furcula presence, pigmentation, post-antennal organ, scales and reproduction type were split into 2 attributes, body length into 10 levels and visual organs (ocelli) into 4 levels ( Table 3).
Because the literature was not exhaustive, we obtained functional traits for only 129 species out of the 143 species, which however represented 98% of the total abundance.

Functional Trait-Based Indices
Regarding the taxonomic approach, these indices were used to assess the structure and the composition of the trait functional communities. The structure of the functional diversity can be assessed through a series of four synthetic indices, which quantify several aspects of functional diversity for a community with species distributed across a multidimensional functional space (Petchey et al., 2004;Villeger, 2008;Mouchet et al., 2010). Functional richness (FRic) represents the functional hyperspace occupied by the community. Functional evenness (FEve) corresponds to the regularity of the distribution of species abundance within this volume. Functional divergence (FDiv) 1 https://portail.betsi.cnrs.fr/   (Laliberté and Legendre, 2010). For the composition of the functional communities, we calculated the community-weighted mean (CWM) for each trait, which is the sum of the taxa affinities to a trait attribute, weighted by the taxon-relative abundance (Lavorel et al., 2008; see formula in Table 4).

Statistical Analysis
To avoid biases induced by pseudo-replicates, we calculated mean species abundances for each site. Even though the number of replicates differs from one land use to another, we assume that the sites are comparable as they were chosen in order to evaluate collembolan communities within these specific land use types. Because of the nested structure of the sampling protocols used in this study, spatial distribution may affect the data distribution. Also, as the sites are distributed around France according to four major climates (oceanic and degraded oceanic, continental, and Mediterranean), the climatic area was considered as a random effect. Generalized linear mixed models (GLMM) were used with a Poisson distribution for the abundances and species richness and functional richness data analyses, as well as a binomial distribution for Pielou's evenness and functional evenness.  Mouchet et al. (2010) and CWM (Lavorel et al., 2008).

Indices Formula
Functional richness N is number of taxa; p i , relative abundance of species i; t ik : trait attribute k of species i; dG: mean distance to the center of gravity; d: sum of abundanceweighted deviances; | d|: absolute abundance-weighted deviances from the center of gravity; PEW: partial weighted evenness; aj is the abundance of species j and zj is the distance of species j from the weighted centroid c. Models were followed by post-hoc tests using the R package "lsmeans" (Lenth, 2016) and lme4 (Bates et al., 2015).
Models were followed by post-hoc tests using the R package "lsmeans" (Lenth, 2016) and lme4 (Bates et al., 2015). For the functional community approach, a matrix was built with CWM values of eight functional traits (129 taxa × 20 CWM calculated on trait categorical levels i.e., attributes of each functional trait). For binary attributes, only one attribute per taxa was used in the matrix in order to avoid auto-correlation. In order to study the effect of land uses on community composition, we performed a redundancy analysis constrained by land use and conditioned by climatic area.
Two partial redundancy analyses (pRDA) were performed in order to evaluate the effects of land use on (i) functional and (ii) taxonomic Collembola communities both with and without climate as a covariable. The statistical significance of the pRDA was assessed by the Monte Carlo Permutation test.

Functional and Taxonomic Indices
The collembolan functional richness (FRic) was high in forest and urban vegetable garden soils (0.19 and 0.21), intermediate in arable land and grassland (0.11-0.16) and low in vineyard and SUITMA soils (0.08, 0.09) ( Table 5 and Supplementary 1). The functional divergence (FDiv) was significantly higher in urban vegetable gardens (0.77) and lower in vineyard (0.57), while other land uses showed intermediate values (0.67-0.65) with no significant differences. The highest value of functional dispersion (FDis) came from forest and urban vegetable garden soils (0.21 and 0.20) and the lowest from vineyard (0.11), the other land uses Frontiers in Ecology and Evolution | www.frontiersin.org

Multivariate Analysis on Functional (CWM) and Taxonomic Composition (Species Distribution)
A Partial Redundancy Analyses (pRDA) was used to test the multivariate effect of land uses on the Collembola taxonomic community and functional trait composition respectively, while climate was used as a co-variable. The CWM of each trait are available in Supplementary 2. Partial RDA independent of climate, showed a significant effect on land uses, on functional composition (13%, p < 0.001) as well on taxonomic composition (13%, p < 0.05). The barycenter of each land use on the first axis of each pRDA is presented in Figure 1 (biplot in Supplementary 3).
We observed on the left of the axis with three groups of land use that discriminate between functional compositions: (i) arable land, (ii) grassland, vineyard soil and urban vegetable gardens, (iii) forest and SUITMA. Alternatively, three groups differed by their taxonomic composition: (i) arable land, grassland and vineyard soil, (ii) SUITMA and forest, and (iii) urban vegetable gardens. However, there is frequent overlap on the pRDA between the two approaches (taxonomic and functional). Only vegetable gardens and the vineyards greatly differ. In contrary, pRDA on taxonomic and functional indices overlap totally for forest. For the other land uses, we observed a part of overlap.

DISCUSSION
The effects of human activities on soil communities usually highlight an anthropization gradient from undisturbed or only slightly impacted soils (duration and intensity of human activities) under forests to more highly impacted urban soils (Joimel et al., 2016). This anthropization gradient has been questioned with regards to soil biodiversity, showing a greater Collembola taxonomic diversity in SUITMAs than in arable soils (Joimel et al., 2017). The present results on collembolan functional diversity, by highlighting a much more complex reality, invite to explore the use of soil functional indices to assess land use effects, in addition to more common indices (taxonomic or geochemical indices).
Indeed, when looking at collembolan functional diversity, SUITMAs appear to be clearly less diversified than when we consider them from a taxonomic point of view. Particularly, we have demonstrated here that the functional richness and evenness of SUITMAs are similar to vineyard soils. So, irrespective of the disturbance factor (urbanization or agronomic stress), their functional diversity was lower than that of forest soil. These results are in agreement with frequently reported expectations that functional diversity indices (richness, evenness, and dispersion) often to decrease along disturbance intensity gradients (Mouillot et al., 2013).
A functional redundancy disturbance hypothesis is suggested in the literature. This states that specialized species are more vulnerable due to filtering of the same traits and that generalist species have greater resistance (Lechêne et al., 2018). This is potentially visible in our results when looking at urban vegetable gardens, where taxonomic diversity is high but the functional composition is similar to that of other land uses.
This biotic homogenization in highly anthropized soils could be related to three different factors. The first of these is soil pollution, due to the increase in metal or organic pollutants in the related sites, which is often highlighted as being responsible reducing the diversity and activity collembolan communities (Fountain and Hopkin, 2004), or inducing shifts in species or functional composition (Cebron et al., 2011;Santorufo et al., 2014;Joimel et al., 2018b).
Secondly, homogenization can change due to an introduction of exotic species (McKinney and Lockwood, 1999), or a selection of synanthropic species, as demonstrated by many previous studies on plants in urban areas as a whole, or in particular types of urban green areas (Marco et al., 2008). As the native status of Collembola species is not clear, our study cannot confirm whether such a functional redundancy is indeed induced by biological invasion. Nevertheless, an increase in the rate of generalist species already present in disturbed sites is highly probable. We noticed an increase in the abundances in urban and industrial land uses of collembolan species described in the literature as "synanthropic, " including Folsomia similis Bagnall, 1939or Coecobrya tenebricosa Folsom, 1902(Joimel et al., 2018aVincent et al., 2018).
Thirdly, homogenization can be impacted by landscape configuration, such as colonization by Collembola in a highly disturbed context like urban green spaces or agricultural landscapes, where they could encounter obstacles due to FIGURE 1 | Graphical display of the first axis of Partial Redundancy Analyses (pRDA) on the functional (13% explained) and taxonomic composition (13% explained) of collembolan community. Functional composition corresponds to community -weighted means (CWM) of 8 functional traits resulting a list of 27 attributes, whereas taxonomic composition corresponds to community assemblage (129 species). Horizontal lines represent the 95% C.I. urban settlement patterns (Joimel et al., 2018b), or be limited by landscape homogenization. Other studies have highlighted the potential importance of the surrounding landscape on Collembola species assemblages .
Although we have taken into account the potential effect of climate, other factors independent of the current land use (previous land use, soil properties influenced by the geochemical context, agricultural practices, and biotic interactions) may also strongly affect the structure of Collembola communities. Further studies are thus needed to differentiate between the effects of these various environmental factors at different scale, in order to better understand global results obtained for each land use.
Conclusions regarding land uses effects will differ depending on the indices used, functional or taxonomic ones. Rigorous and widely applicable indicators of biodiversity are needed to monitor the responses of ecosystems to global change (Vandewalle et al., 2010). The functional approach leads to a new series of bioindicators for soils, which would be complementary to the physico-chemical or taxonomic indices currently used. The difference between urban vegetable gardens and other SUITMA displayed with the trait analysis, shows that the trait-based approach could be used as a complementary tool when survey soil biodiversity, which again could be used alongside the commonly used species-based approach. Even if species richness can reveal a pattern for land use (Joimel et al., 2017), components of functional diversity may provide supplementary information on biodiversity changes thanks to their mechanistic standpoint (Martins da Silva et al., 2015). Functional traits could be advantageous as part of a monitoring program on national or international scales, due to the independence of biogeographic regions. Indeed, on a broader scale, taxonomic richness could be strongly correlated to pedoclimate.
Moreover, the development of a functional trait-based approach, independent of any taxonomic identification, are especially interesting for Collembola. Indeed, a lower level of knowledge could be less required to measure functional traits, than for identifying Collembola at the species level. The development of trait measurements instead of morphological identification could be done only with soft, easy-to-measure traits and following a standard protocol (Moretti et al., 2016). However, the measurement of some functional traits can be more timeconsuming than morphological species identification, especially for species which are easily identified through binoculars. It must be noted that in the present study collembolans were first taxonomically identified and then related to trait-values using information provided in the literature.
It should also be noted that this study was carried out on communities of Collembola extracted with the Macfaydentype apparatus according to the ISO 23611-2 standard. Other extraction methods exist such as pitfall traps or Berlese. All these methods may lead to an under-or over-estimation of certain Collembola species in function of their ecomorphological groups, Frontiers in Ecology and Evolution | www.frontiersin.org 6 July 2021 | Volume 9 | Article 630919 depending in particular on the use or not of light, in addition to a heat gradient.

CONCLUSION
Thanks to our data on 3,056 samples of Collembola communities in 788 sites, we provided a range of values for different taxonomic and functional indices in France and for a broad range of land uses. The complementary results using taxonomic and functional parameters argue in favor of promoting, not only a wide taxonomic diversity, but also a more functional biodiversity, in order to survey soil biodiversity in various land uses.
Our study also illustrates the value of integrating urban and industrial land use studies in advancing our knowledge of soil biodiversity. This improved knowledge of trait and species assembly in soil fauna communities, as well as its dependence on current land use, will therefore contribute to the management and preservation of soil biodiversity in urban and agricultural land.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. This data can be found here: https://doi.org/10.15454/UU2FQT (DATAverse INRAE, available online with publication).

AUTHOR CONTRIBUTIONS
SJ and JC wrote the manuscript. All authors participated in the acquisition of functional trait or/and abundance data and reviewed the article.