Kin competition overrules spatial selection as driver of range expansions

With ongoing global change, life is continuously forced to move to novel areas, which leads to dynamically changing species ranges and imposes rapid shifts in biotic communities and ecosystem functioning. As dispersal is central to range dynamics, factors promoting fast and distant dispersal are key to understanding and predicting species ranges . During range expansions, genetic variation is strongly depleted and genetic homogenisation increases. Such conditions should reduce evolutionary potential, but also impose severe kin competition. Although kin competition drives dispersal, we lack insights into its contribution to range expansions, relative to other processes, such as dispersal evolution. To separate evolutionary dynamics from kin competition, we combined simulation modelling and experimental range expansions using the spider mite Tetranychus urticae. Both modelling and experimental evolution demonstrated that plastic responses to kin structure increased range expansion speed by about 20%, while the effects of evolution and spatial sorting were marginal. This insight resolves an important paradox between the loss of genetic variation and earlier observed evolutionary dynamics facilitating range expansions. Kin competition may thus provide a social rescue mechanism in populations that are forced to keep up with fast climate change. . CC-BY-NC-ND 4.0 International license not peer-reviewed) is the author/funder. It is made available under a The copyright holder for this preprint (which was . http://dx.doi.org/10.1101/150011 doi: bioRxiv preprint first posted online Jun. 14, 2017;


Introduction
The speed and extent of range expansions and biological invasions have traditionally been regarded as consequences of human introductions or ecological factors such as enemy release (Keane & Smith 2002).Recently, a booming field of theory has demonstrated the importance of evolutionary dynamics through spatial selection and sorting of dispersal and traits related to reproduction at the expanding range front (Shine et al. 2011).The process of genetic assortment at expanding range borders results in the evolution of increased dispersal because highly dispersive genotypes colonize vacant habitat first.
In addition, systematically low densities at the leading edge select for increased reproductive performance.Emerging assortative mating from spatial sorting then accelerates these evolutionary dynamics at the range front (i.e. the Olympic village effect (Phillips et al. 2010)), speed up range expansions and make biological invasions even more challenging to contain (Phillips 2015).
Paradoxically, spatial sorting of genotypes during invasion is tightly associated with a successive loss of genetic variation due to subsequent founder effects, rendering drift important, and thus potentially ).Dispersal is recognized as a central and independent trait in life history (Bonte & Dahirel 2017), and known to have a strong impact on spatial dynamics in fragmented landscapes or during range expansions (Kubisch et al. 2014;Legrand et al. 2015;Cheptou et al. 2017).Even when dispersal entails high costs, dispersers may be released from local competition, thereby increasing their inclusive fitness (Hamilton & May 1977).Theory predicts that high dispersal rates are evolving in metapopulations when relatedness within local populations is high (Duputié & Massol 2013), but relatedness is also known to prompt fast plastic changes in dispersal kernels in organisms that show kin-recognition (Bitume et al. 2013).Kin competition may thus be a key driver of fast range expansions and could potentially explain the paradox of fast expansions despite severe genetic diversity loss (Estoup et al. 2016).To date, quantitative insights on its relevance relative to evolutionary dynamics are lacking.
A series of studies has used experimental range expansions to document evolutionary divergence in life history traits between core and edge populations in protists (Fronhofer & Altermatt 2015), beetles (Ochocki et al. 2017; Weiss-Lehman et al. 2017) and plants (Williams et al. 2016b).The latter studies included reshuffling or replacing experiments to quantify the eco-evolutionary loop, i.e. how evolution feedbacks on the range expansion dynamics.By following this experimental procedure, individuals at the front are systematically replaced by individuals from a source population, or from a random patch in the experimental range expansion to destroy any signature of evolution.These experiments have, however, as a major drawback that also genetic structure is destroyed.Given the earlier outlined relevance of kin competition for dispersal, and the central role of dispersal for range expansions, observed differences in expansion speed can equally be attributed to changes in the genetic structure per se, rather than to intrinsic evolutionary dynamics.Kin recognition is after all not restricted to animals, and has been demonstrated in several plant species as well (Dudley & File 2007; Dudley et al.

2013
), including the model species used by (Williams et al. 2016b).
The impact of relatedness on dispersal kernels have been extensively studied in the two-spotted spider mite Tetranychus urticae Koch (Acari, Tetranychidae) (Bitume et al. 2013).This mite species has been developed as a model in experimental evolution (Macke et al. 2011;Fronhofer et al. 2014;De Roissart et al. 2016;Van Petegem et al. 2016).Spider mite life history traits, including dispersal, are documented to be under selection but equally to be highly plastic in response to inter-and intragenerational environmental and social conditions (Magalhães et al. 2009; Bitume et al. 2011Bitume et al. , 2014;;Fronhofer et al. 2014;Van Petegem et al. 2015).The effect of genetic relatedness on both dispersal distance and emigration rate is, for instance, as strong as that of density dependence (Bitume et al. 2014).Following these insights, we used a simple but highly parameterised simulation model and two series of experimental range expansions to test the relative effect of kin competition and evolution (here combined effect of spatial sorting and spatial selection) on range expansion dynamics.The latter was accomplished by contrasting evolved trait divergence and the rate of range expansion in two sets of experimental range expansions that differed in the level of genetic variation and spatial structure.
While a first set of experiments followed the earlier used replacement manipulations that eliminate both kin structure and evolution, a second set of experimental range expansion prevented only evolution while maintaining kin structure along the gradient.

General model algorithm
The model is individual-based and object-oriented and simulates demographic and evolutionary processes along a one-dimensional array of patches.Patches contain resources, which are consumed by individuals at different rates depending on their life stage (juvenile or adult).Resources are refreshed weekly.We parameterised relatedness-dependent dispersal according to earlier research (partly published in (Bitume et al. 2013).Relatedness-related dispersal was only studied for average mite densities, so no further density-dependence was implemented.We thus assume the relatednessdependent kernels to be relevant for the mean population densities during range expansion.A detailed model description and additional results on in silico trait evolution are available in SUPPLEMENTARY MATERIAL 1.
Males and females of Tetranychus urticae differ in a number of aspects.Firstly, males are much smaller when reaching the adult life stage, and hence contribute little to resource consumption (males hardly feed when becoming adult).Secondly, dispersal behaviour differs between the two sexes, with adult females being the dominant dispersers, whereas juveniles and males disperse very little.Lastly, the species is characterized by a haplodiploid life cycle, where non-mated females only produce haploid male offspring, and mated females can produce both haploid male and diploid female eggs.The sex ratio of spider mites is female-biased, with approximately 0.66 males to females.For these reasons and for the sake of simplicity, we designed the model to only include female mites, where the genotype of the individual is passed on from mother to daughter.Individuals carry the following genetic traits: age at maturity, fecundity, longevity, and a categorical neutral genotype (one unique allele per individual) which defines relatedness.Mean relatedness of an individual A in a patch X can be calculated as the number of individuals in patch X carrying the same relatedness genotype as individual A, divided by the total number of individuals, and hence ranges from 0 (no related individuals present) to 1 (all individuals are related to individual A).After 80 steps, concurring with 80 days in our experimental range expansions, both the length of the metapopulation (i.e. the total number of patches in the metapopulation) and the mean life history trait values at the core and edge were recorded.To this end, individuals present in the first patch of the metapopulation (core) or in the last three occupied patches (edge) were tracked (cf. the experimental part of the study) and the mean value of every life history trait was calculated and recorded.
The following scenarios were tested:

T. urticae strains
Several different strains of T. urticae were used within the current study: LS-VL, MR-VP, SR-VP, JPS, LONDON, and MIX.The LS-VL strain was originally collected in 2000 on rose plants in Ghent (Belgium) and since then maintained on common bean (Phaseolus vulgaris, variety Prélude) in a laboratory environment (Van Leeuwen et al. 2004).This strain is known to harbour sufficient genetic diversity for studies of experimental evolution (Van Leeuwen et al. 2008;De Roissart et al. 2016).The MR-VP, SR-VP, JPS, ALBINO and LONDON strains, in contrast, were collected from different greenhouses and completely devoid of any genetic variation by consistently inbreeding mothers with sons over seven generations (see Díaz-Riquelme et al. 2016 for the followed procedure) .These are the non-evolving kin lines abbreviated as ISO.By crossing these five different isofemale strains, we created an evolving kinline containing substantial genetic variation, further abbreviated as MIX.This was done by reciprocally crossing males and females of each of the isofemale strains: for each combination of strains, one female (last moulting stage) of strain X/Y was put together on a bean patch with three males of strain Y/X, allowing fertilisation (in case a fertilisation was unsuccessful, this step was repeated).The resulting F1, F2, and F3 generations were again mixed in such a manner that, eventually, we obtained one mixed strain (MIX) that comprised a mixture of all isofemale strains.Stock populations of the LS-VL and MIX strain were maintained on whole common bean plants in a climate-controlled room (28.1°C ± 2.1°C) with a light regime of 16:8 LD, while stock populations of the ISO strains were maintained on bean leaf rectangles in separate, isolated incubators (28°C, 16:8 LD).Before eventually using the mite strains to initialise the experimental metapopulations, they were first synchronised.For each strain, sixty adult females were collected from the respective stock populations, placed individually on a leaf rectangle of 3.5  4.5 cm 2 , and put in an incubator (30°C, 16:8 LD).The females were subsequently allowed to lay eggs during 24 hours, after which they were removed and their eggs were left to develop.Freshly mated females that has reached the adult stage on day prior to mating of the F1 generation were then used to initialise the experimental metapopulations (see below).As all mites were kept under common conditions during this one generation of synchronisation, direct environmental and maternally-induced environmental effects (Macke et al. 2011) of the stock-conditions were removed.

Experimental range expansion
An experimental range expansion consisted of a linear system of populations: bean leaf squares (2  2 cm 2 ) connected by parafilm bridges [81 cm 2 ], placed on top of moist cotton.A metapopulation was initialised by placing ten freshly mated one-day-old adult females on the first patch (population) of this system.At this point, the metapopulation comprised only four patches.The initial population of ten females was subsequently left to settle, grow, and progressively colonise the next patch(es) in line through ambulatory dispersal.Three times a week, all patches were checked and one/two new patches were added to the system when mites had reached the second-to-last/last patch.Mites were therefore not hindered in their dispersal attempts, allowing a continuous expansion of the range.A regular food supply was secured for all populations by renewing all leaf squares in the metapopulation once every week; all one week old leaf squares were shifted aside, replacing the two-week-old squares that were put there the week before, and in their turn replaced by fresh patches.As the old patches slightly overlapped the new, mites could freely move to these new patches.Mites were left in this experimental metapopulation for approximately ten generations (80 days) during which they could gradually expand their range.

Treatments
We performed two experiments, in each of which experimental metapopulations were each time assigned to one out of two different treatments (see Figure 1).

Fig 1. The four treatments of the two parallel experimental range expansion experiments.
In the first experiment (upper panels), microcosms were either assigned to NMP (non-manipulated) or to RFS (replacement from stock).Both these treatments harboured a relatively high amount of standing genetic variation (mites from LS-VL strain), but in RFS all adult females were regularly replaced (red crosses) by females from the LS-VL stock while this was not the case for NMP.In the second experiment (lower panels), microcosms were either assigned to an evolving kin (MIX) or non-evolving kin line (ISO).The former harboured standing genetic variation (mites from MIX strain; different isofemale lines represented by a single mite colour), but the latter did not (here only one setup with a single isofemale line represented).No reshuffling was performed in this second experiment.By the termination of the experiment, the final range size was measured as the number of occupied patches ((dashed line thus denotes variable number of patches between the fixed core and emerging edge patch).
In the first experiment, we contrasted unmanipulated control LS-VL strains, further abbreviated as NMP, with a treatment where females in the metapopulations were replaced on a weekly basis by randomly chosen, but similarly aged, females from the LS-VL stock.This Replacement From Stock treatment is further abbreviated by RFS.The metapopulations within this treatment thus started with a high enough amount of standing genetic variation for evolution to act on.Kin structure was not manipulated in this treatment and kin competition was therefore expected to increase towards the range edge (see introduction).As a result, any spatial sorting of phenotypes was nullified and kin structure randomized and hence no longer expected to increase towards the range edge.The latter treatment maintains age and population structure (i.e., if x females were on a patch before the replacement, they were replaced by x females from the stock) but prevents genetic sorting, and destroys local relatedness, thus preventing both spatial sorting and kin competition.In this experiment, we thus compared unmanipulated, genetically diverse metapopulations (NMP treatment) with regularly reshuffled metapopulations where only effects of density-dependent dispersal remained (RFS treatment) (cf. (Ochocki et al. 2017)).Both treatments were replicated six times.
In the second experiment, we compared non-evolving with evolving kin lines (MIX versus ISO), an approach earlier followed to avoid evolution in experimental evolution (Turcotte et al. 2011;Ochocki et al. 2017;Wagner et al. 2017).Experimental range expansions starting with evolving kin lines (MIX), harboured standing genetic variation on which evolution could act.No manipulations of kin structure were performed.Experimental metapopulations with non-evolving kin lines (ISO) were initialised using mites from the SR-VP, JPS or LONDON isofemale strain (originally, there were also metapopulations for the MR-VP and ALBINO strain, but these experimental metapopulations collapsed very early within the experiment).These metapopulations therefore harboured no standing genetic variation for evolution to act on.As in the MIX treatment, kin structure was not manipulated.In this second experiment, we thus compared unmanipulated metapopulations (MIX treatment) with metapopulations where only condition dependency (density-dependent dispersal and kin competition) played a role (ISO treatment) [cf.(Wagner et al. 2017)].Both treatments were replicated six times (in case of ISO, two replicates (i.e., experimental range expansions) per isofemale strain were set up).
Genetic trait variation as determined in common garden experiments did not differ among any of the lines, likely due to the dominance of plasticity [see SUPPLEMENTARY MATERIAL 2].Starting from the same levels of trait variation, MIX and NMP thus represent treatments where range expansions are determined by evolution and putative kin interactions, ISO represents the treatment where kin structure is high but where evolution is restricted, and RFS a treatment with both kin competition and evolution constrained.In addition to monitoring range expansions along the linear system, we quantified life history trait variation and genetic variation in gene expression between core and edge populations (details in SUPPLEMENTARY MATERIAL 3, 4)

Results
We used the highly parameterized, but simple simulation model based on spider mite life histories and relatedness-dependent dispersal reaction norms to formalise our hypotheses and predictions.We simulated one-dimensional range expansion during over 8-10 generations (80 days, SUPPLEMENTARY MATERIAL 1).Despite the incorporation of uncertainties regarding condition-dependent dispersal thresholds, the model predicted range expansion to proceed at a 25.9% slower mean rate when signatures of both kin competition and spatial sorting were removed, while expansion rates were only 7.4% slower when spatial sorting was prevented, but kin competition was present (Fig. 2).Thus, range expansion speed increased in our model with 21% by kin competition, but only by 1% due to spatial sorting

Fig 2. Predicted range size (last patch occupied in the linear gradient) by a highly parameterised, stochastic model simulating expansion dynamics in the experimental setup over a period of 80 days (see supplementary material 1). Overall, range expansions proceed slower when kin competition is excluded (B). Median values (indicated by the red lines) under the scenario with kin competition and spatial sorting (A) are similar to those for scenarios with kin competition but where spatial sorting is prevented (C). The scenario 'No kin competition and spatial sorting' could neither be modelled in the individual based models nor experimentally validated.
To test this prediction we ran two parallel experiments where we started experimental range expansions with a limited amount of founders (10 females), thereby mimicking ongoing range expansion of T. urticae along a linear patchy landscape.Each replicated population invaded its respective landscape for ten generations (spanning 80 days).As outlined above, two parallel experiments were conducted in which genetic diversity and relatedness were manipulated to infer the relative importance of spatial sorting and kin competition for range expansion dynamics.

Fig 3: Observed range expansion averaged (±STDEV: coloured belt) over the six replicates per treatment, population densities are shown along the gradient from core to edge (distance). From generation to generation, the populations advance their range (densities along the linear metapopulation).
For each replicate, we measured range expansion dynamics and genotypic trait structure at an unprecedented level of detail.By counting the number of adult females three times per week on each of the occupied patches during the experimental range expansion, we detected a 28% lower rate of range expansion in the treatment with mites replaced from stock, in which kin competition and evolution were constrained, versus the non-manipulated control treatment (NMP-RFS contrast: day × treatment interaction, F1,54.8=7.62;P=0.007; Fig. 3).However, no statistically significant differences were found between the treatments that contrasted the evolving with non-evolving kin lines, which inhibited spatial sorting but left kin competition intact (MIX-ISO contrast: day × treatment interaction, F1,71.1=0.71;P=0.40).Differences in slopes were 0.082 ± 0.026 SE patches/day for the NMP-RFS comparison and 0.030 ± 0.036 SE patches/day for the MIX-ISO comparison, with eventual range size We subsequently tested whether increased range expansion resulted from evolved trait differences between edge and core populations [see SUPPLEMENTARY MATERIAL 3].No significantly higher dispersal rates were detected in individuals from the expanding front, relative to the ones collected from the core patches.Therefore, the accelerated range advance in the treatments with unconstrained evolutionary dynamics was achieved independently of evolutionary changes in dispersiveness.Consistent with predictions of enhanced intrinsic growth rates leading to faster range expansions and sorting of traits at the expansion front (Burton et al. 2010), we found evolved differentiation in the intrinsic growth rate.Intrinsic growth rates were systematically higher in edge relative to core populations in experiments that allowed evolutionary dynamics (NMP: F1,153=5.32,p=0.0225;MIX: F1,235=6.46,p=0.0117;See Fig. 4), but not in those where evolution was experimentally inhibited.

Discussion
The destruction of spatial genetic structure, and thus of both kin competition and evolution, resulted in a lower expansion rate than just the inhibition of spatial sorting by depleted genetic variation.In all treatments except RFS, kin competition was high due to serial founder effects following small population sizes at the initiation of the experiments.Evolution following spatial selection or spatial selection (Shine et al. 2011) was not pronounced in the experimental range expansions where individuals were not replaced or where genetic variation was not depleted.The impact of evolved differences in growth rates on range expansions were therefore only marginal relative to the impact of eliminating the potential for kin competition.The stronger inhibition of range expansion in the replacement treatment (RFS) did not result from lower levels of trait variation and can therefore only be attributed to the elimination of kin ).Surprisingly, we found no indications of variation in any single fitness-related trait among the different treatments.We observed significant replica*location variation in traits during experimental evolution, eventually resulting in divergent trait covariances among replicates (see SUPPLEMENTARY MATERIAL 3).Under such conditions, different life history strategies encompassing multivariate trait correlations but leading to similar high population growth rates might eventually be spatially sorted.Bootstrapping of the vital rates within replicates could however, not confirm the empirically determined higher growth rates at the leading edges (see SUPPLEMENTARY MATERIAL 3).We therefore attribute this opposing evidence from simulated relative to observed intrinsic growth rates to the fact that the simulated ones systematically assume invariant life trait expression during population growth, so neglecting density effects and other individual interactions.Because the interpretation of fitness should be tested under relevant, and varying realistic demographic conditions (Bonte et al. 2014), our empirical assessment thus approaches realistic conditions better than theoretically composed ones.
We here can exclude heterosis as a driving mechanism leading to a "catapult effect" and subsequent faster range expansions in the MIX treatment (Wagner et al. 2017) since metapopulations were initialised with an already mixed strain instead of with separate strains that would hybridise after initialization.We neither did find differences in quantitative genetic (SUPPLEMENTARY MATERIAL 3) and transcriptomic (SUPPLEMENTARY MATERIAL 4) trait variation between the inbred, mixed, and highly diverse stock population, nor between core and edge populations, again indicating the dominance of plasticity over evolution for life history trait expression in our model system.Ambulatory dispersal in mites has a genetic basis (Yano & Takafuji 2002) but is simultaneously highly dependent on differences in density, also those experienced in earlier generations by parents and grandparents (Bitume et al. 2014).We assessed the mites' dispersal behaviour under conditions that reflect the low-density conditions in higher spread rate variance, we did not detect any changes in spread rate variation.Hence, no general conclusions on the impact of spatial sorting on spread predictability can thus be made as they will likely depend on many joint ecological, evolutionary and social factors.
The obliteration of spatial genetic structure, and thus both kin competition and spatial sorting, resulted consequently in a lower expansion rate relative to experiments were only spatial sorting was prohibited by depleting genetic variation.Our results provide the first evidence of kin competition as an overlooked but quantitatively highly significant driver of range expansions compared to spatial sorting.Emerging genetic structure and high relatedness per se along a range expansion front can thus be responsible for fast range expansions, even in the absence of substantial sorting of individual life history traits.The loss of genetic variation during range expansions and biological invasions is typically considered to be a limiting factor because it constraints the potential for local adaptation.We here show that, on the contrary, that it may actually lead to faster range expansions, impose social rescue and therefore allow population to keep pace with high rates of climate change.(2002).Cooperation and competition between relatives.Science (80-.)., 296, 2015): Offspring~max(0, fecundity_intercept-0.45*age_as_adult).

5) Energy loss: individuals use up part of their resource reserves, depending on their current life
stage.Juveniles lose a daily amount of 3.5 from their resource reserves, adults a daily amount of 7.

6) Resource dynamics forcing:
In order for the dynamics to match the empirical experimental evolution (see suppl.material 2), the metapopulation is extended every two days with new patches so that two empty patches are present at the range expansion front.The model keeps track of the entire length of the metapopulation and is built to track the dynamics as experimentally implemented (see suppl.material 2).Replenishment of resources: Weekly renewal of available resources in patch.The total amount of resources in each patch is replenished by 400 at the end of every 7 days.
After initialization of the model, this time step was run 80 times (cf. the approximately 80-day period of the experimental part of the study).
Models were run in two phases: One run to gain general insights into the relative impact of the treatments on changes in range expansion.Here, we randomly assigned values for reserve-dependent dispersal and reproduction from the range 2-5.
One run with optimised reserve-dependent dispersal to fit observed data on range expansion under the control (evolutionary scenario A).The threshold for dispersal and reproduction was found to be 3, which matches empirical observations that dispersal is highly dependent on body condition.A value of 3 corresponds with the availability of reserves to survive three days of starvation.Subsequently, the optimised reserve-dependent dispersal model was run for 10 000 iterations for scenarios A-C to get output on range expansion rate and life history trait evolution.After 80 days, both the length of the metapopulation (i.e. the total number of patches in the metapopulation) and the mean life history trait values at the core and edge were recorded.To this end, individuals present in the first patch of the metapopulation (core) or in the last three occupied patches (edge) were tracked (cf. the experimental part of the study) and the mean value of every life history trait was calculated and recorded.

Range expansion rate
The rate at which range expansion occurred was highest under the scenario A, where both relatedness and trait values where passed on to the offspring (mean expansionscenA=27.37 patches), followed by scenario C, where only relatedness is maintained, but trait values are not (mean expansionscenC=25.16 .

Life history trait evolution
Differences in all recorded life history traits (fecundity, longevity, and age of maturity, see figure 2) showed very little differences under all 3 scenarios.This again corresponds with the experimental part of this study, where we found significant changes in population-level population growth rate, but where distinct evolution of individual level life history traits could not be detected.The most apparent trend in the model output is the difference in spread in trait values between scenario A and B on the one hand, and scenario C on the other hand.Given that scenario C is the only scenario where trait values are consistently reset for all offspring, this seems to indicate that significant drift occurs in scenario A and B. Strikingly, the effect of drift does not appear to differ strongly between the edge and the core, whereas it is usually argued that drift should act much stronger on populations at the edge, due to low population sizes.Results: No differences in CV for trait variation could be detected across traits for each of the used lines (see figure X.1; ANOVA: F4,43=0.71;P=0.59).For the separate traits, consistent variation was not detected either (Table X.1).Overall, isofemale line S was characterised by the highest level of standing phenotypic variation prior to the experimental evolution.Egg survival, juvenile survival, development time, and lifetime fecundity we measured at the individual level for multiple offspring from a single dam and variance could be partitioned among and within dams.Using a full-sib design, we estimated the dam component of the explained variation of lifetime fecundity as 0.15 in the stock-LSVL, 0.45 in the mix, and 0 in all isofemale lines.The dam component of longevity was lower with 0.18 in the stock-LSVL, 0.03 in the MIX, and 0 in all isofemale lines.According to our expectations, phenotypic variation among relative to within dams was thus highly reduced in the isofemale lines.
Figure S3.1:Genetic trait variation averaged across all traits for the different used lines.
Table S3.1:overview of mean values, standard deviation (SD) and coefficient of variation (CV) for each of the measured traits in the different lines.

Range expansion dynamics
Three times a week, the number of patches (populations) within each metapopulation and the number of adult females within all populations of these metapopulations were counted.As such, we kept track of the position of the range front (after ten generations of range expansion, providing us with a measure of the rate of range advance) and attained additional information regarding population dynamics (local densities, hence density curves).
Trait Divergence between core and leading edge Prior to the assessments of trait evolution, mites were first sampled from the core and edge patch of each experimental metapopulation.For each core and edge, ten adult females were collected, placed with five together on a bean leaf rectangle of 3.5  4.5 cm 2 and put in an incubator (30°C, 16:8 LD).The females were subsequently allowed to lay eggs during 24 hours, after which they were removed and their eggs were left to develop.The resulting synchronised freshly mated one-day-old adult females of the F1 generation were used for the trait assessments described below.As all mites were kept under common conditions during this one generation of synchronisation, direct environmental and environmentally-induced maternal effects [cf.(Macke, Magalhaes et al. 2011)] of the experimental metapopulations were levelled out.

Ambulatory dispersal behaviour
The ambulatory dispersal behaviour of the mites was tested using a linear system (cf.experimental metapopulations) of four bean leaf squares (2  2 cm 2 ) connected by parafilm bridges of 1  8 cm 2 .All experiments were run in a climate-controlled room (28.1°C ± 2.1°C; 16:8 LD), after which ten one-day-old adult females were each time put together on the first patch of the linear system.During one week, the position of the ten mites was then recorded daily.Afterwards, emigration/immigration rates were assessed as the number of mites leaving the first patch/entering the last patch over time, and the mean distance dispersed, the maximum reached mite densities on the last patch of the linear system, the (first) day of maximum mite densities on the last patch, and the day of the first mite arrivals on the last patch were calculated.

Intrinsic rate of increase (r)
The intrinsic rate of increase [r, see (Birch 1948)] was assessed by putting a one-day-old adult female on a bean leaf square of 5.5  5.5 cm 2 and allowing her to establish a population.For each edge or core, six Petri dishes (each with one bean leaf square with a single female) were put in an incubator (30°C, 16:8 LD).From day 8 onwards, the number of adult females was counted weekly (for NMP and RFS) or twice a week (for MIX and ISO) during three weeks.
Life-history traits that affect r: egg survival, juvenile survival, development time, sex ratio, adult size, and daily and lifetime fecundity and longevity In a first setup, egg survival, juvenile survival, development time (for females and for males), sex ratio, and adult size (for females only) were assessed following the exact same logic and procedures as described in Van Petegem et al. (2016).However, we assessed sex ratio and adult size within the same setup as egg survival, juvenile survival, and development time.In the previous study, a separate setup was used for sex ratio and for adult size.Therefore, all the traits were assessed for the offspring produced by one fresh adult female during 24 hours of egg-laying and developed at 20°C.In the previous study, sex ratio was measured for all offspring produced during the first 7 days of egg-laying and developed at 27°C, while egg survival, juvenile survival, and development time were assessed starting with three four-day-old adult females.The females were each put on a bean leaf rectangle of 1.5  2.5 cm 2 (which is slightly smaller than in the previous study -only for practical reasons), of which three were put in each of two (for NMP and RFS) or five (for MIX and ISO) Petri dishes.Sex ratios of 1 (i.e., all males, indicating that the mother was not inseminated) were not included in the data analysis.In a second setup, daily fecundity, lifetime fecundity, and longevity were assessed following the exact same logic and procedures as described in Van Petegem et al. ( 2016), but the females were each put on a bean leaf rectangle of 1.5  2.5 cm 2 (cf.above), of which three were put in each of three Petri dishes and stored in an incubator at 30°C [27°C in Van Petegem et al. (2016)].Both mean and cumulative daily fecundity were assessed.For mean daily fecundity, days on which the female became stuck into the cotton (bean patches were always placed on moist cotton) or died were excluded from the data analysis.

Leaf consumption
Leaf consumption (as a measure of foraging) was calculated by comparing photographs of the bean leaf rectangles from the second setup (see above) taken at day 0 (before the female was placed on her patch) with photographs taken at day 3 or day 5. Photographs were taken with a Nikon D3200 camera positioned at a fixed distance from the surface to be photographed (bean leaf).To obtain high-quality photographs (in terms of colours), two LED lights illuminated this surface from a 45° angle through a blue coloured light filter.All photographs were analysed using ImageJ 1.47v (Wayne Rasband, National Institutes of Health, Bethesda, MD, USA), by running a script that enabled us to calculate the surface of the bean leaf that was consumed between day 0 and day 3 or between day 0 and day 5.

Univariate statistical analyses
Most analyses were performed in SAS 9.4 (©2013 SAS Institute Inc., Cary, NC, USA), using either the GLIMMIX (generalised linear mixed model, for binomially distributed data) or MIXED (linear mixed model, for normally distributed data) procedure.In the first case, parameters were estimated via pseudo-likelihood techniques, while in the second case, restricted/residual maximum likelihood (REML) was applied.For all generalised linear mixed models, a binomial error structure was modelled with the proper link function (logit link) and potential overdispersion was corrected for by modelling residual variation as an observation-level random factor (Verbeke and Molenberghs 2000).Effective denominator degrees of freedom for the tests of fixed effects were always computed according to a general Satterthwaite approximation.Because the variance explained by random effects varied among the different dependent variables in our study, these effective denominator degrees of freedom were different for each statistical model.
For the assessment of range expansion dynamics (more specifically the position of the range front), comparisons were made between the treatments in each of our two experiments (NMP vs. RFS and MIX vs. ISO -see above) and treatment was thus set as the independent fixed effect.Data were analysed as repeated measurements, with residuals following a first order autoregressive covariance structure (lowest AIC relative to other covariance error structures).
For all assessments of trait evolution, in contrast, comparisons were made within treatments.Here, the origin of the mites (core vs. edge) was the independent variable.In those cases where we were interested in a linear trend over time (rate of range advance, cumulative daily fecundity, intrinsic rate of increase (r), emigration/immigration rate, and mean distance dispersed), day and its interaction with mite origin were added as two additional independent fixed effects.The analyses were then restricted to that part of the data showing a linear trend (i.e., for cumulative daily fecundity, analyses were performed only for those days when the mother mite still produced eggs; for intrinsic rate of increase, analyses were restricted to the days before reaching carrying capacity of the bean leaf; for emigration/immigration rate and mean distance dispersed, analyses were only run for those days prior to reaching an ideal free distribution of the mites).For the intrinsic rate of increase (exponential curve), the data were thus first log-transformed.All longitudinal data were analysed as repeated measurements, using covariance models with autoregressive covariance error structure (lowest AIC relative to unstructured, Toeplitz, or exponential errors).
For fecundity, and longevity, offspring from single dams were put on single leaf rectangles within petri dishes, thereby enabling us to control for within-dam variation.For survival and development time, up to three offspring from dams were placed on single leaf rectangles, and therefore, we also added leaf rectangle identity (nested within Petri dish) as a random effect to control for dependency among the mites developed on the same leaf.Furthermore, for all assessments of trait evolution, replicate was modelled as a random effect to control for dependency among the mites originating from within the same replicate (there were six replicates (i.e., metapopulations) for the NMP, RFS, and MIX treatment and three times two replicates for the ISO treatment -see above).In case of the ISO treatment, replicate was moreover nested within strain (JPS, SR-VP or LONDON -see above).Accordingly, for all assessments involving mites from the ISO treatment, strain was included as an additional random effect to control for dependency among mites originating from the same ISO strain.
Finally, for all assessments of trait evolution, we also added the interaction between replicate and mite origin as an extra random effect to take into account that the effect of mite origin might differ between replicates (e.g. in replicate 1, mites from the core patch could be more fecund than mites from the edge patch(es), while the opposite trend might be found in replicate 2).
Bootstrap of synthetic life history part Differences between core and edge leaf squares in performance were assessed using the basic reproductive rate R0, corrected for generation time, as a metric following the approach of Cameron et al. (Cameron et al. 2013).Briefly, R0 was estimated from the basic life history data using the following formulas: R = S × P × F is the number of adult females produced per female per generation, with S the survival rate to maturity, P the proportion of females among individuals surviving to maturity and F the lifetime egg production.T = Ti + Tm is the average cohort lifespan sensu Cameron et al. (Cameron et al. 2013), with Ti the age at which maturity is reached and Tm, the average time from an individual's maturity to the birth of one of its offspring.

R0 = exp(ln(R)/T)
As the life-history traits underlying the calculation of R0 were obtained from different individuals, we combined them by bootstrap resampling (10 000 iterations).This enabled us to obtain mean values and confidence intervals for R0 at each replicate × leaf square type (core/ edge) combination.R0 at core and edge leaf squares were then compared for each replicate in each treatment, and were deemed different if 95% confidence intervals around means were non-overlapping.

Ordination of the mean trait values across replicates
A Principal component analysis of trait variation among replicates and position was performed to detect directional responses of trait expression among the microcosm experiments.

Covariance of traits across replicates
We produced heat maps of trait variance/covariance structure based on trait spearman correlation among replicates within edge or core populations of each of the four treatments.Consistent patterns of covariation are informative on putative genetic correlations among traits (i.e., syndromes).Non-consistent covariance structure indicates trait associations emerging within each of the treatments.

Results
Trait Divergence between core and edge populations Ambulatory dispersal behaviour In the first experiment, no differences were found between edge and core mites in terms of immigration rate, emigration rate, mean distance dispersed, maximum mite densities on the last patch of the linear system, the (first) day of maximum mite densities on the last patch, and the day of the first mite arrivals on the last patch for both the NMP and RFS treatment.In the second experiment, mean distance dispersed (F1,92.2=20.90;p<0.0001) and the number of mites immigrating into the last patch of the setup (F1,91.8=22.85;p<0.0001) both increased more rapidly for edge compared to core mites for medium inter-patch distances in the ISO treatment, but not in the MIX treatment (no differences found).No differences in mean distance dispersed and immigration rate were found between edge and core mites for short or long inter-patch distances for both the ISO and MIX treatment.Furthermore, no differences were found between edge and core mites in terms of emigration rates, the maximum mite densities on the last patch of the linear system, the (first) day of maximum mite densities on the last patch, and the day of the first mite arrivals on the last patch for all inter-patch distances in both the MIX and ISO treatment.

Intrinsic rate of increase (r)
In the first experiment, edge populations exhibited a significantly higher intrinsic rate of increase compared to core populations for the NMP treatment (F1,153=5.32;p=0.0225), while no differences were found for the RFS treatment.In the second experiment, edge populations exhibited a significantly higher intrinsic rate of increase compared to core populations for the MIX treatment (F1,235=6.46;p=0.0117), while no differences were found for the ISO treatment.
Life-history traits that affect r: egg survival, juvenile survival, development time, sex ratio, adult size, daily and lifetime fecundity and longevity In the first experiment, the cumulative number of daily produced eggs increased significantly slower in the edge compared to the core for the NMP treatment (F1,635=12.83;p=0.0004), while no differences in cumulative daily fecundity were found for the RFS treatment.Furthermore, egg survival was significantly lower (F1,6.486=5.84;p=0.0490) and adult size significantly higher (F1,131=6.07;p=0.0151) in the edge (egg survival: 0.8998 ±0.0320; adult size: 98432 ±2564.70)compared to the core (egg survival: 0.9778 ±0.0120; adult size: 94992 ±2576.40)for the RFS treatment, while no differences in egg survival and adult size were found for the NMP treatment.No differences in juvenile survival, development time (for both females and males), sex ratio, mean daily fecundity, lifetime fecundity and longevity were found between edge and core of the NMP or RFS treatment.In the second experiment, the cumulative number of daily produced eggs increased significantly faster in the edge compared to the core for both the MIX (F1,847=51.54;p<0.0001) and ISO (F1,1049=5.48;p=0.0194) treatment.Furthermore, juvenile survival was significantly higher in the edge (0.9513 ±0.0179) compared to the core (0.8399 ±0.0381) for the ISO treatment (F1,5.284=10.77;p=0.0202), while no differences in juvenile survival were found for the MIX treatment.No differences in egg survival, development time (for both females and males), sex ratio, adult size, mean daily fecundity, lifetime fecundity and longevity were found between edge and core of the MIX or ISO treatment.

Leaf consumption
In the first experiment, no differences in leaf consumption were found between edge and core of the NMP or RFS treatment.
In the second experiment, leaf consumption after three days was significantly lower in the edge (3.5657 ±0.4306) compared to the core (4.8391 ±0.4575) for the MIX treatment (F1,64=4.11;p=0.0.0469), while no differences in leaf consumption were found after five days, nor after three days for the ISO treatment.
Estimates of trait variation within each of the lines at the start of end position are available on DRYAD [data sheets can be provided to reviewers].

Bootstrap of synthetic life history part
As can be witnessed from the table S4.1 and Figure S4., none of the bootstrapped population growth parameters were significantly different (non-overlapping upper-lower 95% credibility intervals) between core and edge patches across the replicates.
constraining further evolutionary change (Weiss-Lehman et al. 2017).This loss of genetic variation increases local levels of genetic relatedness (Newman & Pilson 1997; Kubisch et al. 2013; Nadell et al. 2016).In many arthropods, for instance, single female colonisers found highly related populations (Dingle 1978).With increasing relatedness, kin-competition will become one of the major interactions (West et al. 2007), thereby impacting population and community dynamics (West et al. 2002; Nadell et al. 2016).One of the best known responses to kin-competition is an increase in dispersal (Bowler & Benton 2005 A. A treatment where dynamics include putative kin competition and evolution.In this scenario, females pass their allele values to the offspring.Mutations occur at a rate of 0.001 and change the trait value to a randomly assigned one as during the initialisation phase.The genotype ID remains unchanged (relatedness is unaffected by trait value mutation).B. A treatment where dynamics do not allow evolution, by changing trait values during reproduction at a mutational rate of 1.In this scenario, all trait values are reset according to the initialisation procedure.Only genotype is maintained, and therefore kin structure and possible kin competition.C. A treatment representing a reshuffling of females.Under this scenario, and as in the experimental procedure, adult females are replaced each week by random females from the stock population.Thus, both trait values and relatedness genotype are re-initialised, eliminating both kin competition and evolution.The model was entirely programmed in Python 3.3.Syntax of the code is publically available on Github: git@github.ugent.be:dbonte/Python-code-Van-Petegem-et-al.-2017.git

Fig. 4 .
Fig. 4. Observed intrinsic growth rate (averaged over 6 replicates) in the two parallel range expansion experiments.The increase in population size over time is shown for populations that were initialised with a female mite collected after 80 days originating from the core versus edge of the experimental range expansions.
competition, and not by different evolutionary trajectories.Such phenotypic variation despite genetic depletion is probably widespread in wild populations and typically maintained by individual differences in development (Cressler et al. 2017), but we can neither exclude long term intergenerational plasticity (Bitume et al. 2014).Evolved increased growth rates at the leading expansion front accord with processes of spatial selection at the expanding front and are in line with recent studies based on field observations (Phillips et al. 2010; Shine et al. 2011; Alex Perkins et al. 2013), including earlier work on natural population of the used model species here (Van Petegem et al. 2016).Systematically low densities at the range front select for increased reproduction, and such strategies are known to advance range expansion substantially (Phillips et al. 2010; Shine et al. 2011; Van Petegem et al. 2016

Figure S1. 1 :
Figure S1.1:Model output of life history trait evolution (as histograms of mean trait values at core and edge).First row = Mean fecundity, Second row = mean longevity and third row = mean age of maturity.Left column = scenario A, middle column = scenario B, right column = scenario C. Blue bins represent population from core, green bins from the edge of the metapopulation.Full blue and green lines represent mean values over all 10000 simulations for core and edge respectively.

Figure S4. 1 :
Figure S4.1:Mean per-replicate R0 as a function of origin (core versus edge) and treatment, estimated by bootstrap.

Figure S4. 2 :
Figure S4.2: ordination plot showing the start-(core, circle) and endpoint (edge, cross) for each replicate of the NMP (green) and RFS (orange) treatment.The blue arrows represent the five most influential traits in explaining the variation along the first and second axis of the PCA.None of these traits, however, showed a clear directional response within the microcosms.There is a high variability in trajectories with arrows pointing in all directions.This could indicate a stochastic effect within replicates (e.g.gene surfing during range expansion within microcosms).. LONG: longevity, DAFE: mean daily fecundity, CUFE: cumulative daily fecundity, LIFE: lifetime fecundity, DIMO: mean distance moved.-1.0 1.0 -0.6 1.0 DIMO

Figure S4. 3 :
Figure S4.3: ordination plot showing the start-(core, circle) and endpoint (edge, cross) for each replicate of the MIX (green) and ISO (orange) treatment.The blue arrows represent the five most influential traits in explaining the variation along the first and second axis of the PCA.None of these traits, however, showed a clear directional response within the microcosms.There is variability in trajectories, though in the MIX treatment, endpoints somewhat cluster together.The starting points are scattered, which could indicate a stochastic effect among replicates (i.e.different founder populations for the different microcosms) (due to the initial founder effect), but see tableS5.LONG: longevity, DAFE: mean daily fecundity, CUFE: cumulative daily fecundity, LIFE: lifetime fecundity, LONG: longevity.

Figure S5. 1
Figure S5.1 Exploratory PCA plots of the core (red) and edge (green) mite populations of the 6 biological replicates within the ECO-EVO treatment.Multivariate data analysis was based on the corrected cyanine intensities (panel A) or on the Box-Cox transformed (using the caret package in R) corrected cyanine intensities (panel B), representingthe transcript levels on a genome-wide scale.PC1 and PC2 represent 43.4% and 18.8% (panel A) or 38.1% and 13.0% (panel B) of the total data variability, respectively.

Table S4 .1: Means and confidence intervals for R0, as estimated by bootstrap
Ordination of the mean trait values across replicates