Bioemco Jerome Mathieu
email: jerome.mathieu@upmc.fr

T Jonathan Davies

Glaciation as an historical filter of below-ground biodiversity

Keywords: Body size, community assembly rules, dispersal, earthworms, latitudinal gradient, past climate, range size, soil biodiversity

Aim

The latitudinal gradient in species richness is one of the most studied biodiversity patterns.

Here we explore a -north-south gradient in earthworm diversity, and evaluate the importance of current and historical filters in shaping the distribution of present-day below-ground species richness.

Location France.

Methods Using high resolution data on earthworm distributions across France, we document the latitudinal gradients in alpha (, beta ( and gamma ( diversity. We relate these gradients to species' traits, taxonomic aggregation and co-occurrence patterns, and correlate them with the present climate and the history of glaciation in Europe.

Results

We found that -diversity decreases from south to north whereas -diversity increases along the same latitudinal gradient. Communities in formerly glaciated regions are composed of smaller, more mobile species and show trait and taxonomical aggregation. In more southerly populations, which did not experience glaciation, earthworm species are larger, have smaller geographical ranges, and communities demonstrate a decrease in species co-occurrence resulting in lower local species richness.

Main conclusions

We show that species richness gradients can present differentsometimes oppositelatitudinal trends depending upon the scale of the analysis. This scale dependence sheds new light on the underlying causes of global biodiversity gradients. The opposing latitudinal trends of the different components of diversity suggest that recolonization following glaciations during the Pleistocene acted as an environmental filter, and that competitive exclusion may be a more dominant ecological force in these former refugial areas. Overall our results show that past climate changes have left a deep footprint on present-day earthworm diversity patterns, from community to macroecological scales, and that different mechanisms of earthworm community assembly may predominate at different latitudes.

INTRODUCTION

At global scales species richness of most major clades peaks in the tropics. The search for mechanistic explanations has typically focused on correlates with the present environment [START_REF] Currie | Energy and large-scale patterns of animal-and plant-species richness[END_REF][START_REF] Francis | A globally consistent richness-climate relationship for angiosperms[END_REF]Hawkins & Porter, 2003a;[START_REF] Buckley | Environmental and historical constraints on global patterns of amphibian richness[END_REF][START_REF] Powney | Hot, dry and different: Australian lizard richness is unlike that of mammals, amphibians and birds[END_REF], calling upon processes such as latitudinal variation in evolutionary rates [START_REF] Rohde | Rapoport's rule is a local phenomenon and cannot explain latitudinal gradients in species diversity[END_REF][START_REF] Allen | Global biodiversity, biochemical kinetics, and the energetic-equivalence rule[END_REF][START_REF] Mittelbach | Evolution and the latitudinal diversity gradient: speciation, extinction and biogeography[END_REF], or ecological explanations, including competition or environmental carrying capacity (see [START_REF] Willig | Latitudinal gradients of biodiversity: Pattern, process, scale, and synthesis[END_REF]. However, mounting evidence suggests that historical processes are also important in structuring biodiversity gradients [START_REF] Jansson | The fate of clades in a world of recurrent climatic change: Milankovitch oscillations and evolution[END_REF]Hawkins & Porter, 2003b;[START_REF] Davies | The influence of past and present climate on the biogeography of modern mammal diversity[END_REF]. Because contemporary environmental variables and historical process, such as the intensity of long-term climate oscillations, covary closely, distinguishing between drivers has proven difficult (e.g. [START_REF] Jansson | Global variation in diversification rates of flowering plants: energy vs. climate change[END_REF]. However, different processes should leave different signatures on the various components of diversity, namely alpha (α), beta (β) and gamma(γ). For example, rapid in situ speciation is predicted to result in high α-diversity plus high β-diversity and/or γ-diversity, whereas, species persistence through adaptive range shifts or migrations might result in more broad-ranged species, maintaining high α-diversity but low β-and/or γ-diversity.

To date, most analyses exploring regional diversity gradients have used range-map data, representing interpolated species distributions from sparse point location data, thereby limiting comparison between the various components of diversity. Recent evidence suggests that environmental correlates of species richness might vary with data type, for example, point location versus gridded datasets [START_REF] Hurlbert | Species richness, hotspots, and the scale dependence of range maps in ecology and conservation[END_REF]. We suggest that these interpolated data might best represent γ-diversity, whilst point location captures α-diversity, perhaps leading to differences in environment-richness relationships between these different components of diversity. Few studies have explored variations in the strength of latitudinal gradients between differing diversity components because of the scarcity of suitable data at both broad and fine spatial scales [START_REF] Meynard | Beyond taxonomic diversity patterns: how do α, β and γ components of bird functional and phylogenetic diversity respond to environmental gradients across France?[END_REF]; but see a recent paper by [START_REF] Kraft | Disentangling the drivers of  diversity along latitudinal and elevational gradients[END_REF].

Here, we decompose the latitudinal gradient for earthworm diversity into its constituent components (α-, β-and γ-diversity). We use a unique point dataset comprising a complete quantitative inventory of all earthworm species in France, representing over 1300 sites evenly distributed across the country (Fig. 1). In addition to its broad spatial extent and high resolution, this dataset is remarkable in its homogeneity in sampling quality as it was compiled in its entirety by M. Bouché [START_REF] Bouché | Lombriciens de France: ecologie et systématique[END_REF] using standardized protocols and taxonomy. Here, we describe the gradient in contemporary earthworm species richness, and relate it to climatic trends over the past 20,000 years.

The glacial history of France is well documented (e.g. [START_REF] Dercourt | Atlas Peri-Tethys of palaeoenvironmental maps[END_REF][START_REF] Buoncristiani | Quaternary glaciations in the French Alps and Jura. Quaternary glaciations -extent and chronology[END_REF][START_REF] Calvet | Recent advances in research on Quaternary glaciations in the Pyrenees. Quaternary glaciations -extent and chronology[END_REF] and the life history of earthworms suggests that current species distributions in previously glaciated regions must be explained by recent (post-glacial) migration and recolonization [START_REF] Bouché | The establishment of earthworm communities. Earthworm ecology: From darwin to vermiculture[END_REF]. Our results therefore provide an example of a biodiversity gradient where historical processes are thought to be important, and we suggest that a better understanding of species' responses to historical climate change might help in predicting future responses, complementing autecology approaches such as niche-based distribution modelling. Understanding earthworm biodiversity in the context of climate changes is important because they play a prominent role in soil functioning and in the maintenance of ecosystems services [START_REF] Bouché | Stratégies lombriciennes. Soil organisms as components of ecosystems[END_REF][START_REF] Lavelle | Earthworm activities and the soil system[END_REF], and, at larger spatial scales, may be a major contributor to global carbon sequestration [START_REF] Bossuyt | Protection of soil carbon by microaggregates within earthworm casts[END_REF].

Based on our understanding of earthworm ecology and the history of climate change in Europe, we make the following predictions.

1. Earthworms will show a traditional latitudinal gradient with fewer species in the north.

2.

Species currently at higher northern latitudes should be good dispersers and hence at these latitudes β-diversity will be low and range sizes large.

3.

Species with northern distributions should be a subset of species found at more southerly latitudes as they have been through an historical filter of glacial history and, as a consequence, earthworm communities in northern latitudes should be more aggregated in their functional traits related to dispersion and range size.

4.

Biotic processes, such as competition, should be more important in structuring communities in the south where there has been a longer history of continuous coexistence.

MATERIALS AND METHODS

Data sets

Species distribution data were obtained from [START_REF] Bouché | Lombriciens de France: ecologie et systématique[END_REF], and comprise the abundance of all earthworm species in France across more than 1300 sites (Fig. 1), using a standardized sampling protocol. Each site represents an area of homogenous land use, typically covering 1 ha or more.

Sites were spaced up to 30 km apart (Fig. 1), and sampling was performed between 1963 and 1968. Each region was sampled in at least two different years. Samples were taken during autumn and spring, the best climatic periods for earthworms. The basic environmental features of all sites, including location, vegetation cover, elevation, soil chemical properties are reported in [START_REF] Bouché | Lombriciens de France: ecologie et systématique[END_REF] Data are available from the Dryad Digital Repository (see Data Accessibility below)).

Earthworm diversity was sampled by collecting three blocks of soil (1 m × 1 m × 0.3 m) at each site, which were sorted with a wet sieving machine designed specifically for this task. Juveniles, which cannot be identified accurately to species, were reared until maturity to allow identification.

We complemented these data by compiling a database of the range maps of French species across Europe, based on the Fauna Europaea database (de Jong, 2013).

For the set of 105 species present in France, we synonymized taxonomy and species names according to the classification of Fauna Europaea (de Jong, 2013) and compiled information on more than 30 anatomical traits, based on [START_REF] Bouché | Lombriciens de France: ecologie et systématique[END_REF] and [START_REF] Sims | Earthworms: keys and notes to the identification and study of the species[END_REF]. We selected traits related to mobility, ecological preferences and reproduction because we a priori expected such traits to be subject to filtering by environment or involved in competition (see Appendix S1 in Supporting Information for the full list of traits. Finally, we obtained data on the limits of permafrost and refugia during the Late Glacial Maximum (LGM) from several sources [START_REF] Dercourt | Atlas Peri-Tethys of palaeoenvironmental maps[END_REF][START_REF] Buoncristiani | Quaternary glaciations in the French Alps and Jura. Quaternary glaciations -extent and chronology[END_REF][START_REF] Calvet | Recent advances in research on Quaternary glaciations in the Pyrenees. Quaternary glaciations -extent and chronology[END_REF].

Climatic data for the present and the LGM were extracted from the ECHAM3 palaeoclimatic model [START_REF] Braconnot | Results of PMIP2 coupled simulations of the Mid-Holocene and Last Glacial Maximum -Part 1: experiments and large-scale features[END_REF] for mean annual temperature (MAT) and mean annual precipitation (MAP). Habitat complexity was characterized as the variance in elevation (Elevvar), measured as the standard deviation in elevation within each region. Elevation data were retrieved from the French National Institute of Geography, at a resolution of 25 m (Bd alti database, http://www.ign.fr/).

Biodiversity gradient analysis

To explore the latitudinal gradient in earthworm biodiversity, we compared α-, β-and γdiversity components. Alpha diversity is the local species richness at each site. Gamma diversity represents the regional species pool, defined as the 150-km radius around each site: a reasonable scale given the geographical distribution of plots and the relatively limited dispersal distance of earthworms [START_REF] Eijsackers | Earthworms as colonizers of natural and cultivated soil environments[END_REF]. To correct for unequal sampling within regions, we used bootstrapping (n = 100 replicates) to randomly select the same number of sites (n = 28 sites) per region, and computed diversity across this subset. Beta diversity (SOR) was estimated within regions following the approach of Balsega and colleagues [START_REF] Baselga | Partitioning the turnover and nestedness components of beta diversity[END_REF], and was decomposed into a spatial turnover component (SIM) and a nestedness component (NES). The relationship between explanatory variables and the different components of diversity was analysed following Baselga [START_REF] Baselga | The relationship between species replacement, dissimilarity derived from nestedness, and nestedness[END_REF], using Pearson's correlations (r) and Dutilleul's correction for the presence of spatial autocorrelation [START_REF] Dutilleul | Modyfying the t-test for assessing the correlation between 2 spatial processes[END_REF].

Species geographical range size

We used the observed distribution of species in Europe as an indicator of species geographical range size. First, the geographical range of each species was calculated by summing the area of the countries in which the species occur, based on distribution maps of Fauna Europaea and the published literature. Second, we determined the average community geographical range size at each sampling point, by computing the average range size of the species present at each location.

Cumulative species plot

We used cumulative species plots to explore the role of past glaciations as filters on present-day communities. First, we calculated the accumulation of species diversity in 1.2° latitudinal bands moving from south to north. Second, we repeated the procedure, but moving from north to south.

If, as we predict, northern species are a subset of species in the south, only a few additional species will be recorded moving from south to north, and the cumulative species plot should initially be steep then shallow or flat. In contrast, there should be an initially much shallower cumulative plot moving from north to south. The further south the intersection between the two curves, the greater the evidence suggesting that the more northern species pool is a subset of the southern pool. In our analysis, the two most southerly bands represent mainly non-glaciated areas (NG), while the more northerly bands were either under discontinuous or continuous permafrost during the LGM.

Variation in sampling effort per band (i.e. number of sites) was corrected by bootstrapping with maximal equal sampling size per band.

Community evolutionary structure

In the absence of a well-resolved phylogeny for earthworms, we characterized the evolutionary structure of communities across and within sites using the species to genus ratio by calculating the mean number of species per genus for each site with more than one species [START_REF] Simberloff | Taxonomic diversity of island biotas[END_REF].

Overall, 105 earthworm species are present in France, distributed across 29 genera and six families. Sites with many species per genus might be considered to be taxonomically or evolutionarily aggregated, whereas sites of equivalent richness but with few species per genus might be considered to be taxonomically or evolutionarily dispersed. In order to test the degree of taxonomical aggregation we used a null model approach which compared the observed ratio of species per genus with expectations from randomly assembled communities with the same species richness.

In addition, we quantified the variation in latitudinal range at the species level that can be explained at different taxonomic levels (see [START_REF] Hof | Phylogenetic signals in the climatic niches of the world's amphibians[END_REF]. We performed variance component analyses (VCA) and analysis of similarity (ANOSIM), with a restricted maximum likelihood approach, to test the significance of the observed pattern using the functions 'lme' and 'varcomp' in the APE package within R [START_REF] Paradis | Analysis of phylogenetics and evolution with R[END_REF]. A large proportion of the species variance explained at higher taxonomic levels would indicate strong phylogenetic structure in latitudinal range.

We evaluated the functional diversity among sites by constructing a distance tree using 32 species characteristics (see Appendix S1). These characteristics included morphological and anatomical traits that are commonly used to identify species, and which are related to functional aspects of earthworm ecology. Anatomical traits mainly included measures concerning the presence and the position of particular organs: it is of interest to note that such traits are believed to have evolved in relation to burrowing and feeding activity, key functional aspects of earthworm ecology [START_REF] Sims | Earthworms: keys and notes to the identification and study of the species[END_REF]. We used a hierarchical ascendant classification (CAH) with Euclidean distance and the Ward algorithm to cluster species. Usefully, branch lengths on the tree represent morphological distances, and can be used to calculate functional divergence between species and trait aggregation within communities [START_REF] Petchey | Functional diversity (FD), species richness and community composition[END_REF].

Competition versus filtering

We explored the functional structure of communities across sites using the distance tree and metrics developed within the ecophylogenetics literature that compare the mean pairwise distance (MPD) among taxa between sites, more typically used to describe the phylogenetic clustering of species [START_REF] Webb | Phylogenies and community ecology[END_REF]. Low MPD suggests under-dispersion (species are more similar); whereas high MPD suggests over-dispersion (species are less similar). We compare the empirical distributions of MPD to a null model generated from randomly shuffling species membership across sites whilst keeping site species richness constant. Over-dispersion is traditionally thought to indicate evidence for competition (similar species displace each other), whereas underdispersion is thought to reflect filtering processes [START_REF] Webb | Phylogenies and community ecology[END_REF]. If, as we predicted, competition is more important in structuring communities in the south, these communities should show a greater tendency towards over-dispersion, whilst communities in the north should be more aggregated because of the historical filter of past glaciations.

Quantification of species co-occurrence patterns

We used the C-score [START_REF] Gotelli | Null model analysis of species co-occurence patterns[END_REF], a quantitative index of species co-occurrence, as an indicator of the strength of competition in communities. The C-score is defined as (Ri -S) × (Rj -S) where Ri and Rj represent the total number of occurrences of species i and j, respectively, and S is the number of shared occurrences. The average C-score, calculated over all unique species pairs, summarizes the pattern of co-occurrence as a single metric. Significance was assessed by constructing random communities (n = 200) at each site, shuffling the species present in the regional pool using the same regions as for the latitudinal diversity gradient analysis above. This approach requires fixing both the number of species by site and the number of occurrence of each species, which is a good compromise between Type I and Type II errors [START_REF] Gotelli | Null model analysis of species co-occurence patterns[END_REF][START_REF] Gotelli | Swap algorithms in null model analysis[END_REF]. We then compared the value of the observed C-score with the distribution of null C-scores to estimate the probability of non-random species co-occurrence.

Geographical representation

As illustration, we generated maps of average community species range, number of species per genus, maximum body size in the community, mean trait dispersion and deviance in C-scores from the null model, interpolating values by punctuated kriging using the cross-validated semivariograms and a weighted linear combination of 15 surrounding data points. Statistical tests were only performed on sampled points.

RESULTS

Latitudinal diversity gradients

Regional earthworm species richness (γ-diversity; Fig. 2a) decreases with latitude (r = -0.59, P = 0.05, Table 1) and increases with LGM temperature (r = 0.56, P = 0.05) and present precipitation (r = 0.35, P = 0.03). Beta diversity, quantified by Sørensen's index, follows trends for γ-diversity (r = -0.75, P = 0.03), with higher turnover at lower latitudes (Fig. 2b), and positive correlation with LGM temperature (r = 0.77, P = 0.01) and present precipitation (r = 0.47, P = 0.01). However, decomposing -diversity into its separate components (Table 1, Appendix S2) reveals that turnover and nestedess demonstrate different trends with environment. Turnover decreases with latitude (r = -0.57, P < 0.01) and present precipitation (r = -0.42, P < 0.01) but increases with LGM temperature (r = 0.65, P < 0.01). In contrast, nestedness increases with latitude (r = 0.31, P < 0.01), present precipitiation (r = 0.3, P < 0.01) and variance in elevation (r = 0.09, P < 0.01), but decreases with LGM temperature (r = -0.42, P < 0.01). Overall, at the regional scale, earthworm species richness is higher, turnover in species composition is greater, and communities are less nested at lower latitudes and where temperatures at the LGM were warmer. However, local species richness (α-diversity; Fig. 2c) shows an unexpected and rarely reported counter-gradient, with diversity increasing towards higher latitudes (r = 0.77, P = 0.02) and decreasing with LGM temperature (r = -0.70, P = 0.01).

Most endemic and narrowly distributed species are found in the south (Appendix S3), while species at higher latitudes have greater latitudinal extents (Fig. 3a, r = 0.39, P = 0.05, from the correlation of range size against latitude, Table 1). In addition, species range is negatively correlated with LGM temperature (r = -0.43, P = 0.02), and present precipitation (r = -0.34, P < 0.01). Interestingly, species range is conserved at the genus level (VCA: 38% of variance, ANOSIM: r = 0.29, P < 0.01, Fig. 3b), but not at the family or order level (VCA = 1 and 14%, ANOSIM: r = 0.004, P = 0.45 and r = 0.43, P = 0.14 for family and order respectively).

Community assemblage patterns

Moving from south to north we find that the cumulative species plot is initially steep, and converges on an asymptote at approximately 46º N (Fig. 4a), indicating that many new species are encountered as we cross southern latitudinal bands while few new species are included as we approach more northerly latitudinal bands. In contrast, the cumulative species plots moving from north to south is initially flat but steepens significantly when reaching non-glaciated latitudes (at approximately 46º N), as additional species not present in the northern sites are picked up (Fig. 4a). The two curves intersect at around 44º N, suggesting that species-rich northern latitude communities above 44-46º N are composed of a subset of species from more southerly latitudes.

By comparing species in the south (below 44° N) to species in the north (above 46° N), controlling for sampling effort, we find that 62% of total species are found only in the south while only 10% are found only in the north. Overall, only 28% of species are present in both the south and the north. Three further lines of evidence provide additional support for this nested relationship. First, although the total number of genera in the north is less than that found in the south, the number of species per genus (taxonomic aggregation) increases towards the north (Fig. 4b,c, Table 1).

Second, body size shows a significant latitudinal gradient (Fig. 5a, Table 1, r = -0.26, P < 0.01, from the regression of body size against latitude), with large-bodied species concentrated in the south [below 44° N: average = 28cm and max = 105cm (several large species); while above 46° N: average = 20 cm and max = 57 cm (only 1 rare large species)], whereas small species (< 20 cm) are ubiquitous south to north. Three, trait dispersionwhich is usually interpreted as an evidence of species filtering by environmental constraintsdecreases from south to north (Fig. 5b, Table 1, r = -0.14, P = 0.05, from the regression of standardized MPD against latitude), thus traits are more aggregated in the north than in the south of France. Further, the correlation with trait dispersion is even stronger for MAT (r = 0.21, P < 0.01) and MATLGM (r = 0.19, P < 0.01).

We therefore suggest that earthworms in the north and where temperatures were colder represent a subset of more southerly species that have been filtered on species traits related to dispersal ability (i.e. body size) and temperature tolerances.

Finally, to evaluate evidence for competition, we used the C-score index of species cooccurrences, which quantifies the degree of species distribution overlap. The larger the C-score, the higher the proportion of potential species pairs found not to co-occur naturally and, by implication, the greater the role of competition in structuring community assembly. Earthworm C-scores increase strongly from north to south (Fig. 5c, Table 1, r = -0.45, P =0.05, for the regression of C-score against latitude), and with temperatures at the LGM (r = 0.48, P = 0.02), suggesting that competition was more important in the south, where communities have been less exposed to climate change and therefore have a longer history of competitive interactions.

Environmental correlations

In general, we found that temperature at the LGM (MATLGM), was as good as, or better than latitude in predicting earthworm diversity gradients (Table 1). In contrast, correlations with present-day temperature (MAT) and environmental heterogeneity (Elevvar), characterized as variance in elevation, were mostly non-significant. Although correlations with present-day precipitation patterns (MAP) were significant, with only a single exception (r = 0.21 and 0.19 for correlations of trait dispersion against MAP and MATLGM, respectively; Table 1), correlation strengths were always higher for MATLGM.

DISCUSSION

Traditional but also unexpected patterns of biodiversity

Regional trends in earthworm diversity match classical diversity gradients with higher diversity towards the tropics. A number of mechanisms have been proposed to explain this remarkably ubiquitous biogeographical trend (see [START_REF] Willig | Latitudinal gradients of biodiversity: Pattern, process, scale, and synthesis[END_REF] for a review). However, by comparing α,  and  components of earthworm diversity we reveal an unusual reverse gradient, with higher α-diversity towards the pole. This reverse latitudinal gradient has been observed in only a few organisms (e.g. [START_REF] Crow | Species-diversity in aquatic angiosperms: latitudinal patterns[END_REF][START_REF] Bolton | Global seaweed diversity -patterns and anomalies[END_REF][START_REF] Skillen | Species richness of the Campopleginae and Ichneumoninae (Hymenoptera : Ichneumonidae) along a latitudinal gradient in eastern North American old-growth forests[END_REF][START_REF] Chown | Hemispheric asymmetries in biodiversity -a serious matter for ecology[END_REF], and has typically been explained by local-scale heterogeneity associated with geography, geology, hydrology, or history [START_REF] Skillen | Species richness of the Campopleginae and Ichneumoninae (Hymenoptera : Ichneumonidae) along a latitudinal gradient in eastern North American old-growth forests[END_REF][START_REF] Willig | Latitudinal gradients of biodiversity: Pattern, process, scale, and synthesis[END_REF]. Our study illustrates a reversal in the latitudinal diversity gradient among the various components of diversity. Previous studies that have simultaneously analysed latitudinal variation in α-, β-and γ-diversity in other organisms [START_REF] Kaufman | The structure of mammalian faunas in the New World: from continents to communities[END_REF][START_REF] Clarke | Spatial patterns of diversity in the sea: bryozoan species richness in the North Atlantic[END_REF][START_REF] Stevens | Geographical ecology at the community level: perspectives on the diversity of New World bats[END_REF] show a generally positive (but sometimes no) association between the three components of diversity. We suggest that past glaciations acting as an historical filter, in conjunction with present climate and topography, explain this unique diversity gradient for earthworms.

Dispersal ability as a filter of biodiversity

Earthworms are unable to persist in permafrost over long periods, such as experienced during the

LGM [START_REF] Holmstrup | The relationship between temperature and cocoon incubation time for some lumbricid earthworm species[END_REF]. Even species that burrow in the ground to avoid frost [START_REF] Nuutinen | Worms from the cold: lumbricid life stages in boreal clay during frost[END_REF] were not able to do so during this period as the deep layers of soil were frozen. In consequence, it is usually accepted that glaciations during the LGM extirpated all earthworms species from northern latitudes [START_REF] Tiunov | Invasion patterns of Lumbricidae into previously earthworm-free areas of northeastern Europe and the western Great Lakes region of North America[END_REF], hence species currently found at these latitudes must have recently recolonized from historical refugia. This hypothesis is consistent with previous interpretations of earthworm distribution, and with the distribution of various other Northern Hemisphere taxa (e.g. [START_REF] Bennett | Quaternary refugia of north European trees[END_REF][START_REF] Hewitt | Post-glacial re-colonization of European biota[END_REF]Hawkins & Porter, 2003b;[START_REF] Petit | Glacial refugia: hotspots but not melting pots of genetic diversity[END_REF][START_REF] Habel | The fourth paradigm pattern of post-glacial range expansion of European terrestrial species: the phylogeography of the Marbled White butterfly (Satyrinae, Lepidoptera)[END_REF][START_REF] Svenning | Could the tree diversity pattern in Europe be generated by postglacial dispersal limitation[END_REF].

Recolonization of northern France required a combination of good dispersal capacity and some degree of niche plasticity. Dispersal ability may therefore have acted as a filter on these species. Theory predicts that filtering results in reduced variability of species traits, referred to as trait aggregation [START_REF] Keddy | Assembly and response rules: two goals for predictive community ecology[END_REF]. As niches and traits are typically phylogenetically conserved, filtering can also result in taxonomic aggregationan increase in the relative number of species per genusbecause species in some genera (possessing beneficial traits) will be favoured over species within other genera (lacking such traits).

Although the phylogeny of earthworms is still not well resolved, and taxonomical issues such as cryptic species may introduce some noise in the data, our results provided four strong lines of evidence suggesting that earthworms were filtered into more northerly, previously glaciated, communities: (1) there is an increase in trait aggregation from south to north;(2) communities show taxonomic aggregation, with a higher ratio of species per genus in the north and many southern species falling into small genera not found within more northern communities;

(3) species in northerly communities are a subset of species in more southerly communities; and

(4) communities in the north are composed of species with wider geographical range size (which demonstrates taxonomic conservatism) than in the south, a pattern which has been reported widely in other taxa [START_REF] Stevens | The latitudinal gradient in geographical range: how so many species coexist in the tropics[END_REF][START_REF] Rohde | Rapoport's rule is a local phenomenon and cannot explain latitudinal gradients in species diversity[END_REF].

Further, we show that diversity gradients correlate most strongly with temperatures at the

LGM, and that correlations with present-day temperatures or environmental heterogeneity were generally weaker or non-significant. Our results indicate strongly that historical climate changes have left a deep footprint on present-day earthworm diversity via selective recolonizations following glacial retreats at the end of the LGM. Interestingly, our results correspond with a recent study on scarab bettles [START_REF] Hortal | Ice age climate, evolutionary constraints and diversity patterns of European dung beetles[END_REF], a group that nest in the ground, and therefore are also sensitive to permafrost. Scarab beetles in the north were also found to be a nested subset of those in the south, and phylogenetically clustered.

Body size as a key trait

Identifying the traits that determine a species ability to colonize or invade new habitats is a challenge. We note that only earthworm species that were able to recolonize the north of France and Europe are invasive in Canada and northern USAregions that were otherwise devoid of indigenous earthworm fauna due to Pleistocene glaciations [START_REF] Hendrix | Earthworms ecology and biogeography[END_REF] suggesting that they possess particular traits that predispose them to range expansion, such as high dispersal capacity.

We show that these species are on average smaller than those species restricted to former refugia.

In earthworms, body size, which ranges from 1.8 to 105 cm in France [START_REF] Bouché | Lombriciens de France: ecologie et systématique[END_REF], is strongly related to demographic parameters [START_REF] Evans | Studies on the relationships between earthworms and soil fertility. IV. On the life cycles of some British Lumbricidae[END_REF][START_REF] Lavelle | Stratégies de reproduction chez les vers de terre[END_REF], an important determinant of species' colonizing capacities at large spatial scales. First, larger species produce larger cocoons deposited deeper in the ground [START_REF] Lavelle | Stratégies de reproduction chez les vers de terre[END_REF], which are less likely to be transported accidentally by other animals or by humans [START_REF] Marinissen | Population dynamics of earthworms in a silt loam soil under conventional and "integrated" arable farming during two years with different weather patterns[END_REF]. Second, body size also differentiates species with respect to their strategy of desiccation resistance during dry months: small species spend this period at the cocoon stage, while large species enter into diapause or quiescence, waiting for autumn to complete their life cycle, which increases their generation time [START_REF] Bouché | Stratégies lombriciennes. Soil organisms as components of ecosystems[END_REF]. Therefore, bigger species require more time to reach maturity and complete their life cycle, and are typically considered to be K-strategists, whilst small species are more r-strategists [START_REF] Bouché | Stratégies lombriciennes. Soil organisms as components of ecosystems[END_REF][START_REF] Satchell | r worms and K worms: a basis for classifying lumbricid earthworm strategies. Soil biology as related to land use practices[END_REF], which is associated with greater invasiveness [START_REF] Ehrlich | Which animal will invade? Ecology of biological invasions of North America and Hawaii[END_REF].

Interestingly, our results run contrary to that predicted by Bergmann's rule, which suggests that physiological constraints lead to larger body sizes in colder climates (e.g. at higher latitudes), although as originally formulated the rule applied to intraspecific variation in endotherms. As suggested by [START_REF] Shelomi | Where are we now? Bergmann's rule sensu lato in insects[END_REF], Bergmann's rule should be applied with caution to ectothermic taxa.

Historical refugia: haunted by the ghost of competition past?

In former refugia, climatic conditions were more suitable for species persistence, explaining higher total regional diversity; however, we find that local () diversity is lower than observed across more northerly sites. In addition, community C-scores reveal a lower index of species cooccurrence, and higher trait dispersion in the south, although trends are relatively weak. These results are consistent with the signal of competition shaping community structure in historical refugia: competition theory predicts a lower rate of species co-occurrence and higher trait variability where competitive interactions are strong. We suggest that the much longer history of competitive interactions in former refugia may have resulted in greater competitive exclusion and a decrease in local species richness in the south. However, other factors, including isolation and diversity of glacial refugia might have also contributed to observed richness patterns in this region. Further work is required to evaluate these hypotheses more fully; for example, more comprehensive phylogenetic information is required to evaluate patterns of co-occurrence among close relatives, and identify cryptic species that are not easily distinguished by morphology.

CONCLUSIONS

Earthworms across France demonstrate two ecological gradients that run counter to classic diversity patterns: α-diversity is higher at more northerly latitudes and body size decreases from south to north. By using data on environment at the LGM, we reveal how consideration of historical process can help in our understanding of present-day diversity patterns. We show that the imprint of glacial history is apparent statistically in gradients of species richness, species range distributions and the aggregation of species traits at both macroecological scales and at the community level. We suggest that past glaciations have acted as an historical climate filter on dispersal ability, resulting in opposing latitudinal gradients for the different components of earthworm diversity. Previously glaciated (more northerly) communities are composed of species with larger geographical ranges, are functionally and taxonomically aggregated, and show low spatial turnover in species composition (low diversity), despite high local species richness (diversity). This historical filtering process is also apparent in the lower regional ( diversity of northern communities, from which poor dispersers (species with larger body sizes) were filtered out. In the north, earthworm communities are assemblages of past invaders, while communities within former LGM refugia may have been structured more by competitive interactions, perhaps over much longer timeframes.

ACKNOWLEDGEMENTS Table 1

Values for Pearson's correlation coefficients (r) between each explanatory variable and the different features of earthworm communities in France.

SOR, intra-regional beta diversity (following Balsega, 2010); SIM, spatial turnover; NES, nestedness; d-C-score, difference between observed and simulated C-score in the null model; Elevvar, variance in elevation; MAPLGM, mean annual precipitation for the Last Glacial Maximum (LGM) period; MATLGM: mean annual temperature for the LGM period; MAP, current annual precipitation; MAT, mean current annual temperature. Bold values are statistically significant (P < 0.05). The associated probabilities following [START_REF] Dutilleul | Modyfying the t-test for assessing the correlation between 2 spatial processes[END_REF] correction for the presence of spatial autocorrelation are shown in brackets. * Body size is significantly related to latitude 2 but not to latitude (r = -0.27, P = 0.11).

FIGURE LEGENDS

Figure 1 Sample sites of earthworms in France and permafrost limits during the Last Glacial Maximum (LGM, c.16 ka, according to [START_REF] Ehlers | Quaternary glaciations -extent and chronologya closer look[END_REF]. Ice, areas continuously frozen and under ice during the LGM; continuous, areas where the soil was continuously frozen during the

LGM; discontinuous, areas where the soil was periodically frozen during the LGM; and refugia, areas that were not glaciated during the LGM.

Figure 2

 2 Figure 2 Decomposition of the earthworm latitudinal diversity gradient within France into (a)

Figure 3

 3 Figure 3 Distribution of earthworm species range sizes across France. (a) Geographical range

Figure 4

 4 Figure 4 (a) Cumulative number of earthworm species moving from north to south (N) and from

Figure 5

 5 Figure 5 (a) Average maximal species body size per earthworm community mapped across

Fig. 2

 2 Fig. 2

Fig. 3

 3 Fig. 3

Fig. 4

 4 Fig. 4

Fig. 5

 5 Fig. 5

We thank Marcel Bouché and Bradford Hawkins for discussion on the results of this work. Much of this work was conducted at NCEAS (National Center for Ecological Analysis and Synthesis, Santa Barbara, California) and at the Biogeography Lab of UCSB (University of California Santa Barbara). This work was funded by the French National Research Agency (ANR) as a part of the project EDISP no. ANR-08-JCJC-0023.

SUPPORTING INFORMATION

Additional Supporting Information may be found in the online version of this article: Appendix S1 List of species' anatomical traits used in this study. Editor: Miguel Araújo

Appendix S2