Miquel Ferrín
email: miquelferrin@gmail.com

Laura Márquez

Henning Petersen

Sandrine Salmon

Jean- François Ponge

Miquel Arnedo

Bridget Emmett

Claus Beier

Inger K Schmidt

Albert Tietema

Paolo De Angelis

Dario Liberati

György Kröel- Dulay

Marc Estiarte

Mireia Bartrons

Josep Peñuelas

Guille Peguero
email: guille.peguero@gmail.com

P De

Trait-mediated responses to aridity and experimental drought by springtail communities across Europe

Keywords: climate change, collembola, functional biogeography, joint species distribution models, shrublands, soil fauna

1. The capacity to forecast the effects of climate change on biodiversity largely rely on identifying traits capturing mechanistic relationships with the environment through standardized field experiments distributed across relevant spatial scales. The effects of short-term experimental manipulations on local communities, however, may overlap with regional climate gradients that have been operating during longer time periods, though to the best of our knowledge, there are no studies simultaneously assessing such long-term macroecological drivers with climate manipulations.

2. We analyzed this hypothesis with springtails (subclass Collembola), a principal group of soil fauna, in a standardized climate manipulation experiment conducted across six European countries encompassing broad climate gradients. We crossed community data (near 20K specimens classified into 102 species) with 22 eco-morphological traits and reconstructed their phylogenetic relationships to track the evolution of adaptations to live at different soil depths, which is key to cope with dryness. We then applied joint species distribution models to investigate the combined effect of the regional aridity gradient with the local experimental treatment (drought and warming) over the assembly of springtail communities and tested for significant trait-environment relationships mediating their community-level responses.

Our results show:

(1) a convergent evolution in all three major collembolan lineages of species adapted to inhabit at different soil strata; (2) a clear signature of aridity selecting traits of more epigeic species at a biogeographic scale, and (3) the association of short-term experimental drought with traits related to more edaphic life-forms.

4. The hemiedaphic condition is the plesiomorphic state for Collembola while the adaptations for an epigeic life would have been secondarily gained. Epigeic springtails are more resistant to drought but also have a higher dispersal capacity that allows them to seek more favorable micro-habitats after experiencing drier conditions. The observed relative edaphization of the springtail communities after short-term experimental drought may thus be a transient community response.

The disparity between macroecological trends and fast community-level responses

after climate manipulations highlights the need of simultaneously assessing long-term and short-term drivers at broad spatial scales to adequately interpret trait-environment relationships and better forecast biodiversity responses to climate change.

Introduction

Global biodiversity patterns are undergoing rapid shifts driven by on-going climate change [START_REF] Blowes | The geography of biodiversity change in marine and terrestrial assemblages[END_REF]. Our capacity to forecast such biodiversity responses rely, however, on detailed knowledge of the processes behind species' population rates and the corresponding community reorganizations. Ecological theory predicts that a continued stress will filter out the species with the most vulnerable combinations of traits leading to changes on the functional structure of communities [START_REF] Mouillot | A functional approach reveals community responses to disturbances[END_REF]. This selection process may be either the result of the new abiotic conditions impairing the population growth of some species (strict environmental filtering) or the outcome of a reduced competitive performance under the new environmental conditions [START_REF] Cadotte | Should environmental filtering be abandoned?[END_REF]. Traits, defined as any morphological, physiological, behavioural or phenological feature measurable at an individual level with impacts on species fitness via their effect on growth, reproduction and survival, are thus posited as the common currency for functional biogeography [START_REF] Mcgill | Rebuilding community ecology from functional traits[END_REF][START_REF] Violle | The emergence and promise of functional biogeography[END_REF]. Because trait selection generally precedes diversity loss, it is therefore fundamental to detect those traits capturing mechanistic relationships with stress factors such as climate change if we are to anticipate shifts on species composition and ultimately on ecosystem functioning [START_REF] Berg | Adapt or disperse: Understanding species persistence in a changing world[END_REF][START_REF] Haddad | Species' traits predict the effects of disturbance and productivity on diversity[END_REF].

Coordinated and distributed multisite field experiments of climate manipulation are among the best tools available to infer these trait-environment relationships at relevant spatial scales under a scenario of rapid climate change [START_REF] Halbritter | The handbook for standardized field and laboratory measurements in terrestrial climate change experiments and observational studies (ClimEx)[END_REF]. The effects of short-term experimental manipulations on local communities, however, overlap with regional climate gradients that have been operating during longer time periods which may lead to divergent ecosystem responses across such gradients that hinder our interpretations [START_REF] Reinsch | Shrubland primary production and soil respiration diverge along European climate gradient[END_REF]. This dependency of observed ecological patterns on multiple processes operating across a hierarchy of space-time scales is a classic issue in functional biogeography that strongly affects our predictive ability because general inferences and extrapolations are too often bogged down in local contingencies [START_REF] Levin | The Problem of Pattern and Scale in Ecology: The Robert H. MacArthur Award Lecture[END_REF][START_REF] Mouquet | REVIEW: Predictive ecology in a changing world[END_REF]. Thus, for the development of the socalled predictive ecology, a key is to incorporate this complexity when assessing the performance of the analytical methods available in the macroecological toolkit [START_REF] Mouquet | REVIEW: Predictive ecology in a changing world[END_REF].

Despite the growing concern over the impact of climate change on biodiversity, there has been a disproportionate focus on aboveground communities disregarding the key role of soil fauna on ecosystem functioning [START_REF] Bardgett | Belowground biodiversity and ecosystem functioning[END_REF][START_REF] Filser | Soil fauna: Key to new carbon models[END_REF][START_REF] Guerra | Tracking, targeting, and conserving soil biodiversity[END_REF]. Among this neglected soil biodiversity there are the springtails (subclass Collembola), an abundant and diverse group of small arthropods whose communities show fast reductions in taxonomic, phylogenetic and functional richness due to drought, and additionally, are also a good proxy of ecosystem functions like litter decomposition [START_REF] Peguero | Fast attrition of springtail communities by experimental drought and richness-decomposition relationships across Europe[END_REF]. Springtails present a diverse array of adaptations related to dwelling in specific soil layers (Figure 1). Smaller collembolan species with unpigmented bodies, shorter appendages, reduced sensory organs and mostly parthenogenetic reproduction are associated with more edaphic habitats [START_REF] Chahartaghi | Sex ratio and mode of reproduction in Collembola of an oak-beech forest[END_REF][START_REF] Salmon | Linking species, traits and habitat characteristics of Collembola at European scale[END_REF], while larger species with pigmented bodies, longer appendages, sexual reproduction and conspicuous mechanic and light receptors such as trichobothria and ocelli respectively, are typical from epigeic environments [START_REF] Salmon | Linking species, traits and habitat characteristics of Collembola at European scale[END_REF]. The species inhabiting the upper soil horizons have a greater drought resistance thanks to a lower cuticular permeability that allow them to resist desiccation [START_REF] Kaersgaard | The importance of cuticular permeability, osmolyte production and body size for the desiccation resistance of nine species of Collembola[END_REF]. On the other hand, living deeper in the soil profile implies a lower exposure to variability in temperature and soil moisture, and additionally, soil-dwelling species also present biochemical, physiological and behavioral mechanisms to cope with drought [START_REF] Holmstrup | Dehydration Tolerance and Water Vapour Absorption in Two Species of Soil-Dwelling Collembola by Accumulation of Sugars and Polyols[END_REF]. Accordingly, a multi-site replicated experiment suggested that increasingly drier conditions may eventually favor edaphic over epigeic springtail species after short-term climate manipulations [START_REF] Petersen | Collembolan communities in shrublands along climatic gradients in Europe and the effect of experimental warming and drought on population density, biomass and diversity[END_REF]. Even though, this relative edaphization may show idiosyncratic responses to short-term climate manipulations since in single-site studies at sub-arctic latitudes epigeic species seem to be favored by experimental warming [START_REF] Bokhorst | Extreme winter warming events more negatively impact small rather than large soil fauna: Shift in community composition explained by traits not taxa[END_REF][START_REF] Makkonen | Traits explain the responses of a sub-arctic Collembola community to climate manipulation[END_REF]. Hence, potentially divergent community responses across large spatial scales may hamper our ability to detect general trends, and in particular to identify those traits behind the current climate-induced community reorganizations. This may be even more difficult if broad climate constraints have had time to differentially shape springtail assemblages across such environmental gradients, although this macroecological dimension, to the best of our knowledge, has never been included in studies assessing soil fauna responses to climate change.

The present study unfolds from a standardized climate manipulation experiment replicated across six natural shrublands encompassing the broad gradients of temperature and precipitation present in Europe [START_REF] Beier | Novel approaches to study climate change effects on terrestrial ecosystems in the field: Drought and passive nighttime warming[END_REF]. Here, we investigated the simultaneous effect of short-term experimental climate treatments (i.e., drought and warming) and the regional long-term climate gradient synthesized in the aridity index, an indicator based on mean annual temperature and precipitation, over the composition and trait distribution of springtail communities. First, we harnessed data from a comprehensive sampling of near 20 thousand springtail specimens sorted out in 102 species [START_REF] Petersen | Collembolan communities in shrublands along climatic gradients in Europe and the effect of experimental warming and drought on population density, biomass and diversity[END_REF] and then we crossed it with 22 eco-morphological springtail traits [START_REF] Salmon | Linking species, traits and habitat characteristics of Collembola at European scale[END_REF] to track across the springtail phylogeny the evolution of the suite of adaptations to inhabit in specific soil layers (Table 1). Finally, we leveraged this community data to build joint species distribution models (JSDM) to test for significant associations between traits with climate manipulation treatments and aridity that could be mediating the assembly of collembolan communities. We specifically addressed the following hypotheses and predictions: (i) the adaptations that allow springtails to thrive at different layers of the soil matrix will show a pattern of repeated convergent evolution across the main lineages of Collembola; (ii) the longterm environmental forcing of the large-scale climatic gradient present across Europe has driven a selection of drought-tolerant springtail species with increasing aridity; (iii) the experimental short-term increase in local dryness and warming will be a major disturbance driving fast responses that eventually favor more edaphic springtail assemblages.

Materials and Methods

Study sites and climatic manipulations

The six shrublands studied comprised most of the European climatic regions (Figure S1 and Table S1). Mean annual precipitation (MAP) at the sites ranged from 544 to 1,263 mm, and mean annual temperature (MAT) ranged from 7.4 to 16.1 ºC. The major types of shrublands present in the study area were: Atlantic heathland (UK, United Kingdom; NL, The Netherlands; DK, Denmark), continental forest steppe (HU, Hungary), and Mediterranean garrigue (SP, Spain; IT, Italy). In each study site, nine 20-m 2 plots were divided into three blocks and randomly assigned one plot at each block to a warming treatment, a drought treatment, or as a control. Warming plots were covered with reflective curtains during the night throughout the year inducing a passive night-time warming, while drought plots were automatically covered by transparent polythene curtains during precipitation events that retracted as soon as rain stopped thus avoiding any warming effect [START_REF] Beier | Novel approaches to study climate change effects on terrestrial ecosystems in the field: Drought and passive nighttime warming[END_REF]. Manipulations of temperature and precipitation were carried out in the same way in all sites, but the timing and duration of the experimental drought were adjusted to the local climate regime (Table S2). The sites were established in 1998 (UK, NL, DK, and SP) and 2001 (HU and IT), and on average the warming treatment produced an increase of 0.5 ºC of MAT and the drought treatment a 35% reduction of soil moisture (Table S2).

Springtail sampling, phylogeny and trait information

Springtails (Subclass Collembola) were sampled sequentially (IT, SP, HU, DK, NL, and UK) from April to July 2003 to equalize mean temperatures across sites. Five quadrats of 1.25 m 2 representative of the most dominant plant species were chosen in each plot. Springtails were sampled from plants and soil surface by suctioning with an adapted vacuum cleaner connected to a fauna trap. Specimens dropped were also recovered by placing small polythene boxes with aqueous benzoic acid beneath the plants. Soil springtails were sampled in 10-cm deep soil cores with a surface of 25 cm 2 taken to the lab and put through high-gradient extraction funnels [START_REF] Gjelstrup | Jordbundens mider og springhaler[END_REF]. These procedures were slightly adapted in UK and IT due to an excess of soil moisture and stones, respectively [START_REF] Petersen | Collembolan communities in shrublands along climatic gradients in Europe and the effect of experimental warming and drought on population density, biomass and diversity[END_REF]. Springtails were identified to species level and some dubious specimens were revised by Drs. L. Dányí, R. Jordana, and E. Mateos.

A total of 19,641 springtail specimens were collected and classified into 102 specieslevel entities (Tables S3) and are conserved in glycerol in H. Petersen's personal collection.

We reconstructed the phylogenetic relationships of springtails by means of sequence data from two genes (cox1 and 28s) gathered from public repositories covering 75% of the species found in our communities (accession numbers in Table S3). We applied Maximum Likelihood and Bayesian procedures to obtain a highly-supported consensus tree that agrees with previous systematic works integrating molecular and morphological information [START_REF] Yu | New insight into the systematics of Tomoceridae (Hexapoda, Collembola) by integrating molecular and morphological evidence[END_REF]. See supplementary methods and figure S2 for further details.

Additionally, we collated data of 22 eco-morphological traits from the ColTrait database [START_REF] Salmon | Linking species, traits and habitat characteristics of Collembola at European scale[END_REF] for 47 up to 64 of the species present in our study sites (Table 1 and Table S4). These traits encompass different dimensions of their ecological niche such as life history (e.g., reproduction strategy), dispersal ability (e.g., leg length) and biotic interactions (e.g., sensory organs and defensive structures). Besides, for each species having trait data available we calculated its Eco-Morphological Index (EMI, [START_REF] Parisi | Microarthropod communities as a tool to assess soil quality and biodiversity: A new approach in Italy[END_REF]. This trait-based index provides a simple continuous metric ranging from 1 to 20 and indicates the level of adaptation of a given springtail species to a specific soil layer. Accordingly, they may be classified as atmobiotic (EMI < 2), epigeic (EMI < 6), hemiedaphic (EMI = [6-8]) or edaphic (EMI > 8, Figure 1).

Data analysis

First, we investigated the evolution of the adaptations that allow springtail species to live in specific soil layers by means of the reconstruction of the ancestral EMI score across the phylogenetic tree. This was done with the function fastAnc of the R package phylosig [START_REF] Revell | phytools: An R package for phylogenetic comparative biology (and other things)[END_REF] and assuming a Brownian-Motion model of evolution of the suite of traits associated to the EMI. We also estimated the correlation between the similarity in EMI and species' evolutionary distance (i.e. the phylogenetic signal) through Pagel's λ. Then, we assessed how the evolution across Collembola of the suite of traits associated with the EMI may lead to specific patterns in the phylogenetic structure of springtail assemblages in response to the drought and warming treatments. Thus, we obtained the mean neighbour taxon distance (MNTD) and the mean pairwise distance (MPD) of each springtail assemblage standardized against 999 community randomizations with the R package picante [START_REF] Kembel | Picante: R tools for integrating phylogenies and ecology[END_REF]. By accounting for both indexes we explored possible non-random patterns occurring at different phylogenetic depths since MNTD can capture changes at the tips of the phylogeny (e.g. selection of sister species sharing specific adaptations) while MPD can detect shifts in the in-depth phylogenetic structure (e.g. selection of species belonging to specific lineages with highly conserved adaptations) [START_REF] Cadotte | Phylogenies in ecology[END_REF].

Second, we assessed the role of aridity and the experimental climate manipulations on the trait-based assembly of springtail communities across Europe. To do so, we carried out a series of generalized linear latent variable models (GLLVM), as implemented in the R package gllvm [START_REF] Niku | gllvm: Fast analysis of multivariate abundance data with generalized linear latent variable models in r[END_REF]. GLLVMs are a type of JSDM that extend the basic generalized linear mixed-effects model to multivariate abundance data incorporating a small number of latent variables accompanied by species-specific factor loadings to model correlations between response variables (i.e., all species in the community matrix). These latent variables have a natural interpretation as ordination axes, but have also the capacity to predict new values controlling for known environmental factors using standard model selection tools [START_REF] Niku | gllvm: Fast analysis of multivariate abundance data with generalized linear latent variable models in r[END_REF].

Additionally, GLLVMs allow us to test for trait-environment associations and the potential effects of biotic interactions. Although the effects of the environment and biotic interactions cannot be teased apart, the inclusion of the latter are a major advance of JSDM compared with previous species distribution tools and 'fourth corner' models that do not include the co-variation of the abundances between all species within a series of sites or communities [START_REF] Poggiato | On the Interpretations of Joint Modeling in Community Ecology[END_REF][START_REF] Warton | So Many Variables: Joint Modeling in Community Ecology[END_REF]. Hence, we first built an unconstrained model (without environmental predictors) with two latent variables and a negative binomial distribution and a log link function. The selection of the best number of latent variables and distribution family followed the recommendations of [START_REF] Niku | gllvm: Fast analysis of multivariate abundance data with generalized linear latent variable models in r[END_REF], which are based on goodness-of-fit and residual diagnostics. Afterwards, we added sequentially the aridity index and the local climate manipulation treatments as environmental constraints and performed likelihood-ratio tests (LRT) to assess the predictive improvement of these increasingly complex models, and if so, we computed the percentage of deviance explained by the best model. The aridity index of each site was calculated by means of a modified version of the inverse of the Gaussen Index (iGI) as: iGI = 1/(MAP/(2*MAT))*100 [START_REF] Reinsch | Shrubland primary production and soil respiration diverge along European climate gradient[END_REF]. This index captures the large-scale gradient of aridity currently present across Europe and that can be traced back to the beginning of the Holocene about 10,000 years ago after the last Younger Dryas cold spell came to an end [START_REF] Hewitt | Post-glacial re-colonization of European biota[END_REF]. Finally, to see whether specific traits were selected by aridity and by the climate manipulation treatments favouring more edaphic or epigeic springtail species, we built a series of GLLVMs including or not the interaction between a given trait and these two environmental variables and tested their significance via LRT [START_REF] Niku | gllvm: Fast analysis of multivariate abundance data with generalized linear latent variable models in r[END_REF].

Results

Comparative analyses of the EMI revealed multiple independent evolutions of the collembolan traits associated to live in specific soil layers (Figure 2). In all the major lineages of subclass Collembola there are groups of species showing either an epigeic habitus (Symphipleona and Entomobryomorpha) or, to a lesser extent, a more edaphic habitus (Poduromorpha), apparently as a result of a convergent adaptive evolution of a distinctive suite of eco-morphological traits as depicted in Figure 1. The high and significant phylogenetic signal of the EMI (Pagel's λ = 0.757) pointed out, however, that once fixed phenotypic reversions in these set of traits within a clade are rare and only in a few cases sister species show contrasting soil layer preferences (e.g Tomocerus minor and Oncopodura crassicornis; TOMMIN and ONCCRAS respectively in Figure 2). The analysis of the phylogenetic structure of springtail assemblages showed that under the experimental drought there was a relative increase in phylogenetic clustering according to the MNTD, i.e. the mean distance of the more closely related species in the assemblage was smaller under drought relative to control conditions (P < 0.01). However, the MPD showed no differences across climatic treatments (P = 0.37), suggesting that under all experimental conditions there were rather similar phylogenetic community structure.

The unconstrained ordination of springtail communities displayed a clear grouping of all springtail assemblages according to the country of origin (Figure 3a). This clustered pattern captured the species turnover among our experimental sites, but additionally, the two latent variables implicitly highlighted the steep aridity gradient that pervades the European geography by grouping together the Spanish assemblages close to those from the Italian and the Hungarian sites. Accordingly, this pattern disappeared when the aridity index was included as a constraining predictor (Figure 3b), thus confirming water availability and mean annual temperature as major drivers of springtail communities across Europe. Indeed, the aridity index alone captured 70.6 % of the whole compositional variability (LRT P < 0.05). Regarding the interaction between traits and the environment, the structuring power of aridity was stronger than that of the short-term effect of the experimental climate manipulations (Figure 4). Among the 22 traits tested for trait-environment associations, 8 showed significant (P < 0.05) relationships with aridity, 6 with the experimental drought and only 3 with the warming treatment based on LRTs between equivalent models with and without the corresponding trait-environment interaction (see Table 2 for a complete list of all significant trait-environment relationships).

At a biogeographic scale, all traits showing significant interactions with aridity indicated a clear trend towards more epigeic springtail assemblages as mean annual temperatures increased and precipitation decreased. Indeed, species with a predominance of sexual reproduction over parthenogenesis, higher numbers of ocelli and longer furcula and antennae were positively selected by increasing aridity (Figure 4). On the other hand, the lower presence and abundance of pigmented species and with photoreceptors (ocelli) in the plots under the experimental drought, along with the increase of species with pseudocelli and with a greater minimum and maximum number of vesicles in this typically edaphic defensive organ pointed out to a relative edaphization of the springtail assemblages after this short-term experimental treatment.

Warming only reduced the presence of species with parthenogenesis and slightly increased the minimum and maximum number of ocelli in the springtail assemblages.

Discussion

The capacity to forecast the effects of climate change on biodiversity largely rely on identifying those traits capturing mechanistic relationships with the environment through standardized field experiments distributed across macro-ecologically relevant spatial scales [START_REF] Halbritter | The handbook for standardized field and laboratory measurements in terrestrial climate change experiments and observational studies (ClimEx)[END_REF][START_REF] Mouquet | REVIEW: Predictive ecology in a changing world[END_REF]. Our results demonstrate, however, that the effects of local climate manipulations on communities of springtails overlap with the regional gradients that have been operating during longer time periods.

The broad aridity gradient present across Europe promoted assemblages of species with more epigeic traits as water availability decreases and temperature increases. Despite this long-term and strong aridity driver, the effect of the drought treatment led to a relative edaphization of springtail communities. We discuss below how this disparity of community responses against short versus long-term climate restrictions result from the trait differences between epigeic and more edaphic springtails and why this also warns out that our interpretations of the trait-environment relationships behind the current biodiversity reorganizations of soil fauna may change or even be reversed as climate change continues.

Evolutionary convergence of edaphic traits in Collembola

Springtail species adapted to inhabit at specific soil strata typically share a common suite of eco-morphological traits and, to some extent, they even share a similar trophic niche irrespective of their phylogenetic affiliation [START_REF] Ponge | Vertical distribution of Collembola (Hexapoda) and their food resources in organic horizons of beech forests[END_REF][START_REF] Potapov | Connecting taxonomy and ecology: Trophic niches of collembolans as related to taxonomic identity and life forms[END_REF][START_REF] Salmon | Linking species, traits and habitat characteristics of Collembola at European scale[END_REF]. Our results complement these earlier findings demonstrating that this is the result of an adaptive convergence that has repeatedly happened in the lineages of the three major Collembolan orders. The ancestral ecology of springtails has been of particular interest in light of the still unknown terrestrialization path followed by the subphylum Hexapoda [START_REF] Ghilarov | L'importance du sol dans l'origine et l'évolution des Insectes[END_REF][START_REF] Van Straalen | Evolutionary terrestrialization scenarios for soil invertebrates[END_REF]. According to our dataset, the hemiedaphic condition would be the plesiomorphic state for Collembola while the adaptations for a truly epigeic life would have been secondarily gained, particularly by Entomobryomorpha and Symphypleona. This conclusion agrees with previous morphological assessments [START_REF] D'haese | Morphological appraisal of Collembola phylogeny with special emphasis on Poduromorpha and a test of the aquatic origin hypothesis[END_REF], thus bridging together the cladistics approach with our phylogenetic analysis.

The development of an epigeic habitus is associated with larger and pigmented bodies with longer appendages, sexual reproduction and more developed sensory organs, while the specialized edaphic life-form typically imply a smaller unpigmented body, shorter appendages, reduced or absent sensory organs along with the development of the ability to reproduce through parthenogenesis [START_REF] Chahartaghi | Sex ratio and mode of reproduction in Collembola of an oak-beech forest[END_REF][START_REF] Salmon | Linking species, traits and habitat characteristics of Collembola at European scale[END_REF].

Despite the repeated convergent evolution of these traits at a broad phylogenetic scale, the transitions between epigeic and edaphic life-forms, however, are rarely reversed within a clade so that sister species usually share a similar habitus and hence soil layer preferences. A previous study at a smaller spatial scale has shown that evolutionary close species tend to co-occur due to their akin ecological preferences and competitive exclusion of sister species may be rare notwithstanding [START_REF] Ponge | Spatial and taxonomic correlates of species and species trait assemblages in soil invertebrate communities[END_REF]. This suggests that the relative phylogenetic clustering observed under the experimental drought according to the MNTD metric may mirror springtail assemblages where sister species are coexisting. The apparent discrepancy with the MPD, which points out to a lack of phylogenetic clustering due to the climate manipulations, may actually arise by the fact that in all major springtail lineages there are species relatively specialized to either an epigeic or edaphic lifestyle. Thus, the presence of species from all evolutionary lineages under all experimental treatments results in similar in-depth phylogenetic community structure and may preclude any phylogenetic imbalance in the resulting assemblages. Therefore, this adaptive convergence may explain why in Collembola phylogenetic diversity measures are less sensitive to species losses or community shifts due to warming and drought than functional richness or other traitbased diversity metrics [START_REF] Peguero | Fast attrition of springtail communities by experimental drought and richness-decomposition relationships across Europe[END_REF].

Long-term climate gradients versus short-term climate manipulations

The current climate gradients present in Europe date back to around 10000 years BP and by 6000 year BP the vegetation patterns already resembled that of today [START_REF] Hewitt | Post-glacial re-colonization of European biota[END_REF].

Here we demonstrate that the increase in aridity across Europe nowadays is coupled with springtail assemblages where species with typically epigeic traits are better represented. These positively selected traits at a macroecological scale include longer furculas and antennas, well-developed light and mechanic sensory organs, and a predominance of sexual reproduction over parthenogenesis, all of them being characteristic adaptations to live in the soil surface or its uppermost horizons [START_REF] Salmon | Linking species, traits and habitat characteristics of Collembola at European scale[END_REF]. These atmobiotic and epigeic species have a higher resistance to drought and thermal stress than specialized soil-dwelling springtails thanks to biochemical and physiological mechanisms that include a lower cuticular permeability, the production of sugars and polyols to regulate internal osmolality and a greater plasticity in the fatty acid composition of their cellular membranes [START_REF] Dooremalen | Acclimation responses to temperature vary with vertical stratification: Implications for vulnerability of soildwelling species to extreme temperature events[END_REF][START_REF] Kaersgaard | The importance of cuticular permeability, osmolyte production and body size for the desiccation resistance of nine species of Collembola[END_REF]. But additionally, their larger body size, longer appendages and fully functional visual organs provide these epigeic species with a greater dispersal capacity [START_REF] Ojala | Dispersal of microarthropods in forest soil[END_REF], which allows them to migrate to avoid adverse conditions and choose more favorable micro-habitats within their home range [START_REF] Chauvat | Foraging patterns of soil springtails are impacted by food resources[END_REF][START_REF] Ponge | Decreased biodiversity in soil springtail communities: The importance of dispersal and landuse history in heterogeneous landscapes[END_REF], ultimately making them more resilient against local disturbances [START_REF] Lindberg | Recovery of forest soil fauna diversity and composition after repeated summer droughts[END_REF][START_REF] Malmström | Life-history traits predict recovery patterns in Collembola species after fire: A 10 year study[END_REF].

In contrast to epigeic springtails, the species adapted to live deeper in the soil profile typically migrate downwards to escape from dryness [START_REF] Detsis | Vertical distribution of Collembola in deciduous forests under mediterranean climatic conditions[END_REF][START_REF] Hopkin | Biology of the springtails:(Insecta: Collembola)[END_REF]. This behavioral difference may be behind the observed relative edaphization of springtail assemblages under the experimental drought. Atmobiotic and epigeic species may have dispersed seeking for better patches nearby [START_REF] Chauvat | Foraging patterns of soil springtails are impacted by food resources[END_REF][START_REF] Ojala | Dispersal of microarthropods in forest soil[END_REF][START_REF] Ponge | Decreased biodiversity in soil springtail communities: The importance of dispersal and landuse history in heterogeneous landscapes[END_REF] after experiencing between 2 up to 4 years of drought, while edaphic species were still there likely retreated below and perhaps migrating across the soil profile tracking the daily variation in moisture [START_REF] Detsis | Vertical distribution of Collembola in deciduous forests under mediterranean climatic conditions[END_REF][START_REF] Hopkin | Biology of the springtails:(Insecta: Collembola)[END_REF]. Indeed, the soil matrix provides a remarkable buffering capacity against environmental variation. However, the edaphic species are generally more vulnerable to desiccation and thermal stress than epigeic collembolans [START_REF] Dooremalen | Acclimation responses to temperature vary with vertical stratification: Implications for vulnerability of soildwelling species to extreme temperature events[END_REF][START_REF] Kaersgaard | The importance of cuticular permeability, osmolyte production and body size for the desiccation resistance of nine species of Collembola[END_REF] in spite of their biochemical and physiological adaptations [START_REF] Holmstrup | Dehydration Tolerance and Water Vapour Absorption in Two Species of Soil-Dwelling Collembola by Accumulation of Sugars and Polyols[END_REF]. The effectiveness of this vertical migration strategy may therefore be limited if the drought episode lasts too long. The response to the warming treatments did not lead to major changes in the functional structure of the communities. Our warming treatment was rather mild just increasing from 0.2 to 0.9 ºC the MAT of the experimental plots (Table S2). Additionally, previous research suggests that the diversity and composition of springtail communities may be fairly resistant and resilient to moderate warming [START_REF] Alatalo | Collembola at three alpine subarctic sites resistant to twenty years of experimental warming[END_REF][START_REF] Holmstrup | Soil microarthropods are only weakly impacted after 13 years of repeated drought treatment in wet and dry heathland soils[END_REF][START_REF] Holmstrup | Functional diversity of Collembola is reduced in soils subjected to short-term, but not long-term, geothermal warming[END_REF][START_REF] Peguero | Fast attrition of springtail communities by experimental drought and richness-decomposition relationships across Europe[END_REF][START_REF] Petersen | Collembolan communities in shrublands along climatic gradients in Europe and the effect of experimental warming and drought on population density, biomass and diversity[END_REF].

Conclusions

There is a solid consensus around the fact that the responses of species to withstand climate change can be grouped around two major strategies: to disperse or to adapt [START_REF] Berg | Adapt or disperse: Understanding species persistence in a changing world[END_REF][START_REF] Jump | Running to stand still: Adaptation and the response of plants to rapid climate change[END_REF]. Springtails, and most likely all soil fauna, are no exception and when experiencing increasingly adverse environmental conditions may either 'move or change' [START_REF] Ponge | Move or change, an eco-evolutionary dilemma: The case of Collembola[END_REF]. To move implies from local dispersion up to distribution range shifts. The effectiveness of this strategy relies, however, first on the buffering capacity of the local microhabitats, and secondly on the mismatch (or not) between compositional changes of communities tracking environmental suitability and the velocity of climate change itself [START_REF] Devictor | Differences in the climatic debts of birds and butterflies at a continental scale[END_REF]. Dispersal estimates for collembolan species are scarce but lie in the range of few (~5) centimeters per week during the favorable season, so that springtail assemblages, at a community level, are expected to select habitats by active movement within circles up to 200 meters [START_REF] Chauvat | Foraging patterns of soil springtails are impacted by food resources[END_REF][START_REF] Ojala | Dispersal of microarthropods in forest soil[END_REF][START_REF] Ponge | Spatial and taxonomic correlates of species and species trait assemblages in soil invertebrate communities[END_REF][START_REF] Ponge | Move or change, an eco-evolutionary dilemma: The case of Collembola[END_REF]. Taking into account that birds and butterflies, which both have notably high dispersal capacities, are already experiencing climatic debts of 212 and 135 km respectively [START_REF] Devictor | Differences in the climatic debts of birds and butterflies at a continental scale[END_REF], we may have serious doubts about the ability of Collembola to keep up with climate change through their active dispersal. Like other organisms such as plants that are unlikely to migrate fast enough to track the rapidly changing climate, adaptation must play an increasingly important role [START_REF] Jump | Running to stand still: Adaptation and the response of plants to rapid climate change[END_REF][START_REF] Ponge | Move or change, an eco-evolutionary dilemma: The case of Collembola[END_REF]. There are hopeful examples of soil fauna showing rapid evolutionary changes in response to climate change [START_REF] Bataillon | A replicated climate change field experiment reveals rapid evolutionary response in an ecologically important soil invertebrate[END_REF], and although the only example with a springtail species reported some degree of ontogenetic plasticity (i.e. steeper reaction norms of developmental rates of juveniles than those of adults), there was almost no sign of local adaptation to geothermal warming [START_REF] Kutcherov | Temperature responses in a subarctic springtail from two geothermally warmed habitats[END_REF]. If so, this calls for urgent studies assessing the adaptive potential of springtails, and more generally tracking population dynamics and functional trait shifts at a community level [START_REF] Bardgett | Belowground biodiversity and ecosystem functioning[END_REF][START_REF] Berg | Adapt or disperse: Understanding species persistence in a changing world[END_REF][START_REF] Guerra | Tracking, targeting, and conserving soil biodiversity[END_REF], specially taking into account that for conservation purposes collembolan functional diversity indices may outperform phylogenetic metrics and better correlate with ecosystem functioning [START_REF] Peguero | Fast attrition of springtail communities by experimental drought and richness-decomposition relationships across Europe[END_REF].

Finally, our study also warns out that our interpretations of the trait-environment relationships behind the current biodiversity reorganizations of soil fauna may change or even be reversed as climate change continues. As pointed out by previous studies [START_REF] Alatalo | Collembola at three alpine subarctic sites resistant to twenty years of experimental warming[END_REF][START_REF] Holmstrup | Soil microarthropods are only weakly impacted after 13 years of repeated drought treatment in wet and dry heathland soils[END_REF][START_REF] Holmstrup | Functional diversity of Collembola is reduced in soils subjected to short-term, but not long-term, geothermal warming[END_REF], our work also suggests that the observed shifts of springtail communities to climate manipulations may be transient in time, with epigeic species showing faster responses probably due to their higher vagility. Thus, the resultant relative edaphization we have observed in the springtail communities of our experimentally drought plots could change and even be reversed after some time, in light of the clear and strong selection towards more epigeic assemblages that the aridity gradient has exerted at a larger spatiotemporal scale. To the best of our knowledge there are no studies simultaneously assessing long-term macroecological drivers with short-term climate manipulations at relevant spatial scales and this study demonstrates how important this is if we are to adequately forecast soil fauna responses to climate change. See the Supplementary Methods and Table S3 in the Supporting Information for further details about the phylogeny and for the species abbreviations. Each trait-environment interaction was modelled with a separate GLLVM maximizing the number of species with trait data available (47 to 64 species). See Table 2 in Drought and warming effects are averages during the study period.

Soil moisture reduction applies to the end of the experimental drought period.

TABLE S3

List of the 102 springtail species identified. See [START_REF] Petersen | Collembolan communities in shrublands along climatic gradients in Europe and the effect of experimental warming and drought on population density, biomass and diversity[END_REF] for further information. Accession numbers for the molecular markers (28S rDNA and cox1 mtDNA) used to infer their evolutionary relationships and the availability of trait data are also shown.

FIGURE CAPTIONS Figure 1 .

 CAPTIONS1 FIGURE CAPTIONS Figure 1. Schematic representation of the morphological and functional adaptations to living in different soil layers by the main phylogenetic lineages of Collembola (Order Entomobryomoprha, Symphypleona and Poduromorpha) along with the corresponding range for the Eco-Morphological Index (EMI). Based on Parisi et al. (2005), Potapov et al. (2016) and Salmon et al. (2014).

Figure 2 .

 2 Figure 2. Phylogenetic reconstruction of the Eco-Morphological Index (EMI) across the clade Collembola. Red to yellow indicate lower EMI values associated to atmobiotic and epigeic species, while green to blue denote higher EMI values related to hemiedaphic and edaphic species. Drawings are placed along with the corresponding species or its closest relative to illustrate their morphological and functional adaptations.

Figure 3 .

 3 Figure 3. Ordination based on a generalized linear latent variable model fitted to the collembolan community data (102 species) without environmental predictors (a) and after controlling for the effects of the aridity index (b). Symbol shape correspond to the countries included in the study (DK, Denmark; HU, Hungary; IT, Italy; NL, Netherlands; SP, Spain and UK, United Kingdom), and colour refer to the experimental climate manipulation applied to each plot. Taxa with the largest factor loadings are included as indicator species. See Table S3 in the Supporting Information for species abbreviations. Inset in panel b shows the result of a likelihood-ratio test between the unconstrained and the constrained model along with the increase in variance explained by the aridity index.

Figure 4 .

 4 Figure 4. Level plot for the empirical significant interactions between collembolan traits and short-term experimental drought and the long-term large-scale aridity gradient present in Europe according to generalized linear latent variable models (GLLVM).

  the Supporting Information for a complete list of interactions and model outputs. Statistical significance was assessed by means of likelihood-ratio tests between equivalent models with or without the trait-environment interaction. Coloured squares show significant interactions (P < 0.05) with (•) denoting marginal significance (P < 0.1). The colour scale offers an indication of the sign and magnitude of the selection exerted over each trait in the collembolan communities.

  FIGURE 1

Table 2 .

 2 Statistically significant interactions between springtail traits and environmental conditions.

	617
	Trait	Drought	Warming	Aridity	P-value	Nº of species
	Parthenogenesis	-0.63±0.67	-1.19±0.66 *	-2.44±1.20 *	0.0087
	Sexual reproduction	0.05±0.43	0.53±0.43 *	0.87±0.81 *	0.0087
	Furcula length (1)	0.44±1.53	-1.53±1.72	5.76±3.43 *	0.0000
	Furcula length (2)	-0.75±1.47	-0.71±1.44	5.02±2.85 *	0.0000
	Furcula length (4)	-0.05±0.69	0.44±0.66	7.78±1.77 *	0.0000
	Antennal length (2)	-0.47±3.45	-0.26±3.79	7.21±6.10 *	0.0002
	Antennal length (3)	-0.42±0.56	0.26±0.50	3.07±1.34 *	0.0002
	Antenna-head ratio (min.)	0.01±0.38	0.24±0.34	2.05±0.76 *	0.0000
	Antenna-head ratio (max.)	0.07±0.35	0.40±0.32	2.12±0.71 *	0.0000
	Ocelli	-0.65±0.62 *	0.58±0.61	0.61±1.57	0.0199
	Ocelli number (min.)	-0.16±0.23	0.36±0.27 *	3.18±0.86 *	0.0018
	Ocelli number (max.)	-0.09±0.28	0.44±0.27 *	1.73±0.75 *	0.0001
	Post Antennal Organ (PAO)	0.27±0.51	-0.31±0.52	-2.54±1.35 *	0.0024
	PAO number (min.)	0.53±0.54	0.26±0.44	1.30±0.51 *	0.0000
	PAO number (max.)	0.46±0.41 *	0.20±0.36	1.15±0.44 *	0.0002
	Trichobothria	0.07±0.48	0.50±0.45 •	1.46±1.19 •	0.0828
	Pigmentation	-0.60±0.59 *	0.53±0.58	-0.09±1.35	0.0055
	Pseudocelli	1.28±0.67 *	-0.31±0.93	-0.82±3.67	0.0031
	Pseudocelli number (min.)	0.53±0.28 *	-0.01±0.24	-0.04±0.58	0.0002
	Pseudocelli number (max.)	0.53±0.28 *	-0.01±0.24	0.01±0.55	0.0002
	Each line corresponds to a fourth corner model based on Generalized Linear Latent Variable Models. Effect
	estimates are followed by its 95% Confidence Interval. Significant and marginal interactions are represented by
	"

*" and "•", respectively. P-values show the result of a LRT against an equivalent model without the corresponding trait-environment interaction. Last column indicates the number of species with data available for each trait and thus included in the models.

Table S2 . Experimental manipulations at the study sites.

 S2

	Site code			UK	NL	DK	HU	SP	IT
	Start of the experiment		1998	1998	1998	2001	1998	2001
	(pre-treatment year)
	First treatment year		1999	1999	1999	2002	1999	2002
	Drought
	Timing			May-September April-July	May-July	May-June	May-June	April-October
								October-November
	Precipitation excluded	25	19	18	22	49	16
	(% of yearly total)
	Reduction	in	soil	45	43	41	23	28	27
	moisture (% of control,
	0-20 cm)
	Warming
	Timing			Year-round	Year-round Year-round	Year-round	Year-round	Year-round
	Increase in MAT (C)		0.2	0.3	0.9	0.4	0.6	0.4

DK, Denmark; HU, Hungary; IT, Italy; MAT, mean annual temperature; NL, Netherlands; SP, Spain; UK, United Kingdom.

Acknowledgements

All authors thank the efforts of the many people who helped with the fieldwork during the VULCAN-INCREASE projects. The present study received the financial support

from the European Research Council through the Synergy grant ERC-SyG-2013-610028 IMBALANCE-P, the Spanish Government grant PID2019-110521GB-I00, the Fundación Ramón Areces grant ELEMENTAL-CLIMATE, and the Catalan Government grant SGR 2017-1005. M.F is supported by an FPI grant from the Spanish Government.

The raw data supporting the results of this study will be available at a public repository upon the acceptance of the manuscript.

2 Note that furcula as presence/absence and its length are also protective features related with defensive evasion. Trait data collated from COLTRAIT database [START_REF] Salmon | Linking species, traits and habitat characteristics of Collembola at European scale[END_REF]. Some continuous or discrete numerical traits (e.g. Length, Ocelli, PAO number, Pseudoccelli) have minimum and maximum values recorded separately. See Table S4 for further information on the traits included in this study.

Growing season

April-September

April-October

April-September

April-September

January-May

October-December

January-May

October-December

Dominant species

Calluna vulgaris

Calluna vulgaris

Calluna vulgaris

Deschampsia flexuosa

Populus alba

Festuca vaginata

Erica multiflora

Globularia alypum

Cistus monspeliensis

Helichrysum italicum

Dorycnium pentaphyllum

MAT, mean annual temperature; MAP, mean annual precipitation.

MATs and MAPs apply to the study period (see Table 2).

Species with relative cover above 10% in the control plots during the study period are listed as dominant species. Supplementary Methods Phylogenetic data. The phylogeny of springtails was constructed with sequence data of two genes (cox1 and 28s), obtained from public repositories (GenBank and BOLD). We found information for 77 species (75.5% of all species recorded), although in some cases we had to use sequences for the closest taxa within genus or family level (Table S3). The two genes were aligned separately with CLUSTAL W (Larkin et al., 2007) was preferred due to the large number of incomplete sequences. The Birth Death model was selected as tree prior. Because the goal of the analyses was to obtain an ultrametric tree, we assigned an arbitrary rate of 1 to the 'ucld.mean' parameter of the relaxed clock. The accompanying programs LOGCOMBINER and TREEANNOTATOR were used to remove the burn-in generations, to combine the results of three independent chains, and to select the optimal distribution of posterior premasters and tree values.

Despite the inherent topological uncertainty of any phylogenetic reconstruction, the tree we obtained was consistent with previous systematic works integrating molecular and morphological information, even in correctly placing the Family Tomoceridae closer to Neelipleona and Poduromorpha than to Entomobriyoidae [START_REF] Yu | New insight into the systematics of Tomoceridae (Hexapoda, Collembola) by integrating molecular and morphological evidence[END_REF] (Figure S2).

FIGURE S2

Ultrametric Bayesian consensus tree of the springtails sampled obtained with BEAST under the assumption of lognormal relaxed molecular clock and substitution rates arbitrarily set to 1. Outgroups were removed before the analysis and root was assigned based on previous analyses conducted with maximum likelihood and Bayesian inference including one bristletail and one dipluran as outgroups. Numbers at the nodes are the posterior probabilities, which are correlated with branch color with red denoting the maximum support for a node. See Supplementary Table S3 for interpretation of the species abbreviations.