Katrien Van Petegem

Felix Moerman

Maxime Dahirel

Emanuel A Fronhofer

Martijn L Vandegehuchte

Thomas Van Leeuwen

Nicky Wybouw

Robby Stoks

Dries Bonte
email: dries.bonte@ugent.be

K Van Petegem

E A Moerman

Fronhofer

Kin competition accelerates experimental range expansion in an arthropod herbivore

Keywords: invasions, spatial sorting, experimental evolution, relatedness, Tetranychus urticae

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Cooperation and conflict are central to understanding organismal interactions and their impact on population and community dynamics [3,6]. In group-living species, related individuals may profit from collective behaviour, but competition among kin may eventually outweigh the potential benefits [7]. Dispersal provides a prominent means to avoid competition with kin and conflict more generally. Even when dispersal entails high costs, dispersers may be released from local competition, thereby increasing their inclusive fitness [5]. This effect is even stronger when dispersal is not a fixed genetic trait, but a conditional response [8,9], for instance conditional on kin structure where individuals plastically adapt their dispersal strategy to current levels of relatedness [10]. Dispersal has been recognized as a central and independent trait in life history, known to have a strong impact on spatial dynamics in fragmented landscapes or during range expansions [9,[11][12][13]. Nonetheless, we surprisingly lack knowledge about the consequences of the interaction between kin structure and conditional dispersal for ecological patterns at large spatial scales, such as range expansions.

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 [14]. Recently, a booming field of theory has demonstrated the importance of evolutionary dynamics through spatial selection of dispersal and/or reproduction at the expanding range front [15]. 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, systematic 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 [16]), speed up range expansions and make biological invasions even more challenging to contain [17]. Although mechanisms behind spatial selection and spatial sorting are different, we refer here to both as spatial sorting for the ease of communication.

Paradoxically, evolutionary change should be slow during range expansions as genetic variation already gets depleted early-on due to subsequent founder effects, rendering drift important, and thus potentially constraining further evolutionary change [18]. Simultaneously, high levels of local genetic relatedness emerge due to reduced population sizes [2][3][4]. In many arthropods, for instance, single female colonisers found highly related populations [19]. Kin competition in combination with an appropriate conditional dispersal response may thus be a key driver of fast range expansions and could potentially explain the paradox of fast expansions despite severe genetic diversity loss [20]. Unfortunately, conditional dispersal related to kin interactions have to date been neglected in the context of range expansions and biological invasions (but see one study on the evolution of unconditional dispersal and kin structure during range expansions [4]).

We tested the relative effect of kin competition and spatial sorting on range expansion dynamics by means of in silico simulations and experimental range expansions, using the two-spotted spider mite Tetranychus urticae Koch (Acari, Tetranychidae) as a model organism. This mite species allows us to assess quantitative life history traits in detail [21][22][23][24]. Spider mite life history traits, including dispersal, are documented to be heritable but highly plastic in response to inter-and intra-generational environmental and social conditions [25][26][27]. The effect of genetic relatedness on both dispersal distance and emigration rate is, for instance, as strong as that of density dependence [10].

Material and Methods

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. Individuals start out with a limited reserve of resources, which they will need to replenish in order to survive. 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:

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

Experimental range expansions

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 [40]. This strain is known to harbour sufficient genetic diversity for studies of experimental evolution [23,41]. 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) . By crossing these five different isofemale strains, a strain containing substantial genetic variation was created: the MIX strain. 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 [24] 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 -cfr. [10]], 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 t o diffe e t t eat e ts. I the fi st e pe i e t, the e e assig ed to eithe NMP o RF" . I the NMP t eat e t efe i g to oa ipulated populatio , e pe i e tal etapopulatio s e e i itialised using mites from the LS-VL strain. 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). In the RFS treatment sta di g fo Repla e e t F o "to k , e pe i e tal etapopulatio s e e also i itialised ith ites f o the LS-VL strain, but all adult females in the metapopulations were replaced on a weekly basis by randomly chosen, but similarly aged, females from the LS-VL stock. 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 spatial structure (local densities) of the metapopulations within this treatment was however maintained (i.e., if x females were on a patch before the replacement, they were replaced by x females from the stock). In this first 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. [29]). Both treatments were replicated six times.

I the se o d e pe i e t, e pe i e tal etapopulatio s e e assig ed to eithe MIX o I"O . I the MIX t eat e t fo i tu e of i ed li es -see above), experimental metapopulations were initialised using mites from the MIX strain. This strain harboured standing genetic variation on which evolution could act. No a ipulatio s of ki st u tu e e e pe fo ed. I the I"O t eat e t fo I"Ofe ale li e , e pe i e tal metapopulations 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. [30]]. Both treatments were replicated six times (in case of ISO, two replicates (i.e., experimental range expansions) per isofemale strain were set up).

In addition to monitoring range expansions along the linear system, we quantified life history trait variation genetic variation in gene expression between core and edge populations (details in SUPPLEMENTARY MATERIAL 3, 4)

Results

First, we formalized our hypotheses by means of a highly parameterized, but simple simulation model based on spider mite life histories and relatedness-dependent dispersal reaction norms. 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. 1). Thus, range expansion speed increased in our model with 21% by kin competition, but only by 1% due to spatial sorting

Fig 1. 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 s e ario No ki o petitio a d spatial sorti g ould either e odelled in the individual based models nor experimentally validated.

Second, 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). 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. To specifically test for evolution during range expansion, we determined quantitative genetic differences in life-history traits between core and edge populations and measured patterns of gene expression.

Fig. 2. 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 (reshuffling). 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 fro the LS-VL sto k while this was ot the ase for NMP. I the se o d experi e t lower pa els , microcosms were either assigned to MIX or to 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 a first experiment, we compared the dynamics of range expansion and relevant life history traits between six replicated non-manipulated experimental range expansions (NMP) and six replicates where single adult females were randomly replaced by similar-aged females from the same stock population (further abbreviated as RFS -Repla e e t F o "to k . The latter treatment maintains age and population structure but prevents genetic sorting, and destroys local relatedness, thus preventing both spatial sorting and kin competition. In a second parallel experiment, six replicated experimental range expansions with either depleted or substantial standing genetic variation were contrasted [28][29][30]. Single, different isofemale lines were used for the experiments with strongly depleted genetic variation (further abbreviated as ISO), and recombined isofemale lines were used for the genetically enriched experimental expansions (further abbreviated as MIX). The enriched populations maintain density, genetic, and phenotypic structure. In the genetically depleted lines, spatial sorting is impossible but relatedness patterns are left intact. Quantitative 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 spatial sorting and putative kin interactions, ISO represents the treatment where kin structure is high but where spatial sorting is restricted, and RFS a treatment with both kin competition and spatial sorting constrained.

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 RFS scenario, in which kin competition and spatial sorting were constrained, versus the unconstrained NMP scenario (day × treatment interaction, F1,54.8=7.62; P=0.007; Fig. 3). However, no statistically significant differences were found between the unconstrained MIX treatment and the ISO treatment, which inhibited spatial sorting but left kin competition intact (day × treatment interaction, F1,71.1=0.71; P=0.40). Differences in slopes were respectively 0.082 ± 0.026 SE and 0.030 ± 0.036 SE patches/day, with eventual range size matching the average model predictions. In contrast to other studies finding evolution leading to either increased [18,29] or reduced variation among replicates [30,31]

distance). From generation to generation, the populations advance their range (densities along the linear metapopulation).

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 [32], 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 spatial sorting, 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. Spatial selection was not pronounced in the experimental range expansions that allowed evolution. 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 RFS treatment did not result from lower levels of trait variation and can therefore only be attributed to the elimination of kin 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 [33], but we can neither exclude long term intergenerational plasticity [26]. 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 [15,16,37,38].

Systematically low densities at the range front select for increased reproduction, and such strategies are known to advance range expansion substantially [15,16,22,39]. 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 simulated one 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 [39], our empirical assessment thus approaches realistic conditions better than theoretically composed ones.

We here can exclude heterosis as a driving mechanism leadi g to a atapult effe t and subsequent faster range expansions in the MIX treatment [32] 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, or 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 [40] but is simultaneously highly dependent on differences in density, also those experienced in earlier generations by parents and grandparents [26] . We assessed the ites' dispe sal behaviour under conditions that reflect the low-density conditions in Bitume's o k [26] for which grandparental environmental conditions did not affect dispersal behaviour. The accelerated range expansion does not result from elevated densities at the front, and thus density dependency in dispersal, neither did we find evolution of increased competitive abilities through enhanced foraging efficiency or increased long-distance dispersal at the range front (SUPPLEMENTARY MATERIAL 3). With exception of the RFS treatment, kin competition is expected to be high due to serial founder effects following small population sizes at the initiation of all experiments. Kin competition causes more [41,42] and further dispersal [10], and is hence predicted to speed-up range expansion. Kin recognition is not restricted to animals, and has been demonstrated in several plant species as well [34][35][36][37][38]. Evidence for kin recognition is specifically accumulating in Arabidopsis [45], the model system used for experimental range expansions by Williams et al. [48]). They only performed a reshuffling experiment and were not able to control for potential local kin interactions. Because kin recognition is widespread and not restricted to animals, is may be a highly important but neglected driver of range expansions for a wide variety of life forms. I o t ast to Willia s et al's o k he e so ti g a o ed a ia e i spread rate and parallel but independent experiments [18,31] that reports higher spread rate variance, we did not detect any changes in spread rate variation. 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.

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.

Detailed description of the results

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). Slowest range expansion happened under scenario B where both relatedness and trait values were reset (mean expansionscenB = 19.76). Compared to scenario C, range expansion was accelerated by approx. 38.5% for scenario A and approx. 27.3% for scenario B. Thus, this output suggests that the relative importance of relatedness (and associated kin avoidance) can impact range expansion more strongly than spatial selection/sorting of life history traits (rel. effect relatedness = accele a i ce C accele a i ce A = . % .5 % = .7 ).

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.

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.

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. 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. ISOsee 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. [50]. 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. [50], 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.

  no significant differences in the variation in spread rate were present among any of the treatments (coefficients of variation in the reached distance with 95% CI based on bootstrapping: NMP: 0.246 [0.147-0.276]; RFS: 0.2424 [0.133-0.279]; MIX: 0.279 [0.207-0.314]; INBRED: 0.220 [0.198-0.248]).

Fig 3 :

 3 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).

Fig. 4 .

 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.

Figure S3. 1 :

 1 Figure S3.1: Genetic trait variation averaged across all traits for the different used lines.

Figure S4. 1 :

 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. 3 :

 23 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 S5. 1

 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 S3 .

 S3

		Stock LSVL	125620.5 11207.21	0.089215
	Mean distance	MIX	1.5		1.322876	0.881917
		ISO_J	0.402778 0.21606	0.536426
		ISO_L	1.5		0.707107	0.471405
		ISO_S	1.033333 0.895669	0.866776
		Stock LSVL	1.146032 0.853232	0.681274
	Egg survival	MIX	0.6889	0.457761	0.66448
		ISO_J	0.9452	0.408383	0.43206
		ISO_L	0.8966	0.500515	0.558237
		ISO_S	0.8077	0.684012	0.846863
		Stock LSVL	0.8807	0.4534	0.514818
	Juvenile survival	MIX	0.8387	0.361826	0.431412
		ISO_J	0.942		0.40472	0.429639
		ISO_L	0.8846	0.520396	0.588284
		ISO_S	0.8413	0.6624	0.787353
		Stock LSVL	0.9533	0.410599	0.430713
	Female	MIX	18.0714	2.404182	0.133038
	development	ISO_J	20.0889	3.05207	0.151928
	(days)	ISO_L	17.9744	3.279309	0.182443
		ISO_S	18.3333	4.468021	0.243711
		Stock LSVL	16.422	2.0766	0.126452
	Male	MIX	16.4167	1.187352	0.072326
	development	ISO_J	18.2		3.511726	0.192952
	(days)	ISO_L	17.1333	2.867332	0.167354
		ISO_S	17.5		2.776672	0.158667
		Stock LSVL	15.7053	2.653	0.168924
	Lifetime	MIX	39.1667	26.791	0.684025
	fecundity	ISO_J	43.8333	37.8885	0.864377
		ISO_L	43.5		37.8885	0.871
		ISO_S	41.4		41.5045	1.002524
		Stock LSVL	41.2167	34.58428	0.839084
	Intrinsic growth	MIX	44.33333 25.59427	0.577314
	rate	ISO_J	35.25		10.27538	0.2915
	0.9	ISO_L	56.25		25.34265	0.450536
		ISO_S	32.66667 32.25937	0.987532
	0.6 Sex ratio Trait variation (CV) (proportion	Stock LSVL MIX ISO_J	35.42857 25.38654 0.3333 0.33 0.3077 0.357526	0.716555 0.990099 1.161931
	males) 0.3	ISO_L ISO_S	0.4167 0.4	0.397494 0.517148	0.95391 1.292871
		Stock LSVL	0.3909	0.45	1.15119
	0.0
	Iso_J	Iso_L	Iso_S	Mix	Stock_LSVL
			Line			1: overview of mean values,
	standard deviation (SD) and coefficient of variation (CV) for each of the measured traits in the
	different lines.
	Trait	Line	Mean	SD	CV
	Adult female size	MIX	117228.8 9063.099	0.077311
	(contour)	ISO_J	97173.38 11826.55	0.121706
		ISO_L	102027.6 7826.355	0.076708
		ISO_S	94045.09 14011.69	0.148989

Table S4 .1: Means and confidence intervals for R0, as estimated by bootstrap

 S4 Ordination of the mean trait values across replicates

	treatment	replicate	origin	R0_mean	R0_lower 95% CI	R0_upper 95% CI
	NMP	1	core	1.189	1.145	1.227
	NMP	1	edge	1.219	1.184	1.254
	NMP	2	core	1.174	1.131	1.207
	NMP	2	edge	1.168	1.136	1.198
	NMP	3	core	1.059	0	1.151
	NMP	3	edge	1.128	1.084	1.163
	NMP	4	core	1.124	1.072	1.166
	NMP	4	edge	1.116	1.074	1.149
	NMP	5	core	1.170	1.127	1.206
	NMP	5	edge	1.173	1.145	1.200
	NMP	6	core	1.192	1.147	1.232
	NMP	6	edge	1.207	1.147	1.254
	RFS	1	core	1.164	1.122	1.198
	RFS	1	edge	1.167	1.130	1.195
	RFS	2	core	1.155	1.118	1.187
	RFS	2	edge	1.175	1.145	1.202
	RFS	3	core	1.110	1.052	1.155
	RFS	3	edge	1.160	1.113	1.201
	RFS	4	core	1.187	1.160	1.210
	RFS	4	edge	1.099	0	1.189
	RFS	5	core	1.192	1.168	1.213
	RFS	5	edge	1.112	1.039	1.174
	RFS	6	core	1.146	1.097	1.187
	RFS	6	edge	1.189	1.148	1.223
	MIX	1	core	1.153	1.122	1.181
	MIX	1	edge	1.153	1.123	1.179
	MIX	2	core	1.187	1.165	1.210
	MIX	2	edge	1.198	1.178	1.218
	MIX	3	core	1.100	1.009	1.146
	MIX	3	edge	1.172	1.144	1.198
	MIX	4	core	1.171	1.144	1.197
	MIX	4	edge	1.146	1.104	1.177
	MIX	5	core	1.204	1.186	1.221
	MIX	5	edge	1.201	1.187	1.215
	MIX	6	core	1.189	1.163	1.214
	MIX	6	edge	1.161	1.138	1.183
	ISO-JPS	1	core	1.155	1.137	1.170
	ISO-JPS	1	edge	1.174	1.158	1.189
	ISO-JPS	2	core	1.158	1.121	1.189
	ISO-JPS	2	edge	1.172	1.139	1.201
	ISO-LONDON	1	core	1.125	1.083	1.160
	ISO-LONDON	1	edge	1.143	1.100	1.181
	ISO-LONDON	2	core	1.146	1.104	1.184
	ISO-LONDON	2	edge	1.151	1.114	1.187
	ISO-SRVP	1	core	1.097	1.039	1.141
	ISO-SRVP	1	edge	1.107	1.059	1.145
	ISO-SRVP	2	core	1.119	1.077	1.152
	ISO-SRVP	2	edge	1.122	1.087	1.153

Acknowledgements

This project was funded by the Fund for Scientific Research-Flanders (FWO; project G.0610.11 and G.018017N). D.B. and R.S. were supported by the FWO research network EVENET; RS by the KU Leuven Excellence Center Financing PF/2010/07. Services used in this work were provided by the Flemish Supercomputer Center (VSC), funded by Ghent University, the Hercules Foundation, and the Flemish government (Department of Economy, Science, and Innovation).

Appendices

Supplementary material 1: Details of the simulation model Supplementary material 2: Quantitative trait variation at the onset of the experimental evolution Supplementary material 3: Contrasting quantitative trait (co)variation between core and edge populations Supplementary material 4: Methods and results transcriptomics Supplementary material 1: details of the model Model parameterization and time step At initialization, a metapopulation is created with 4 patches, each with 400 resources available. Within the first patch, 10 founding female adults are created with trait values for age of maturity, fecundity intercept, and longevity. Trait values are assigned based on earlier experimental work [22,23] and this study. Life history trait values are sampled from Gaussian distributions at initialisation, and passed on to offspring conditional on a mutation threshold. The following trait values were assigned:  Age at maturity ~N(7.5, 0.98);  Fecundity_intercept ~ N(7, 0.4);  Longevity ~N (13,1).

The basic time step in the model corresponds with a single day in the experiment and consists of a series of operations executed in the following chronological order: 1) Dispersal: adult females have the possibility to disperse one or two patches in the setup, or to remain in their current patch. Dispersal is dependent on relatedness according to earlier data reported in [10], pilot experiments in preparation of the model, and stochastically implemented as the probability of moving either one (P1) or two (P2) patches (no dispersal mortality is implemented), where: Supplementary material 2: Quantitative trait variation at the onset of the experimental evolution Methods: We collected mites, in the same way as for the the experimental metapopulations, from the inbred (ISO F, ISO J, ISO S), mixed (MIX), and stock lines (STOCK LSVL). Mites were placed in groups of five on a bean leaf square 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 (inbred) mated one-day-old adult females of the F1 generation were used for the quantification of the following traits for single females: mean ambulatory dispersal distance, intrinsic rate of increase (r), egg survival, juvenile survival, development time, sex ratio, adult size, and lifetime fecundity. Since all mites were reared under common conditions during this one generation of synchronisation, direct environmental and environmentally-induced maternal effects are removed. For longevity and fecundity, dam components of the overall variation could be estimated through mixed modelling. We quantified coefficients of variation (CV) for all traits and interpreted this as the level of standing phenotypic trait variation.

Supplementary material 3: Contrasting quantitative trait (co)variation between core and edge populations

General methodology

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. [START_REF] Macke | Experimental Evolution of Reduced Sex Ratio Adjustment Under Local Mate Competition[END_REF])] 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 [START_REF] Birch | The intrinsic rate of natural increase of an insect population[END_REF]] 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.

Supplementary material 4: Methodology and results transcriptomics

Prior to transcriptome analysis, mite populations were sampled from the core and edge patches of the NMP and MIX treatment and reared under common environmental conditions. For the NMP and MIX experiments, we collected six and five independent biological replicates, respectively. RNA samples were extracted from 100-150 pooled 1-to-2-day-adult female mites using the RNeasy mini kit (Qiagen) and were subsequently DNase treated (Turbo DNA-free, Ambion). After running an aliquot on a 1% agarose gel, the quality and quantity of the RNA were finally assessed by a TapeStation instrument and corresponding software (version A.01.05 (SR1)) (Agilent Technologies). All samples had a RIN of 7 or higher. RNA was labelled with cyanine dyes, as described in [START_REF] Dermauw | A link between host plant adaptation and pesticide resistance in the polyphagous spider mite Tetranychus urticae[END_REF]. RNA of the edge and core samples was labelled with cy5 and cy3, respectively. Hybridization of cRNA samples was performed as described in [START_REF] Dermauw | A link between host plant adaptation and pesticide resistance in the polyphagous spider mite Tetranychus urticae[END_REF]. On every array, a cy3-labelled core sample was mixed with the cy5-labelled edge sample of each biological replicate. Slides were scanned with an Agilent Microarray High-Resolution Scanner and raw data were extracted with Agilent Feature Extraction software using the GE2_107_Sep09 protocol (available at the Gene Expression Omnibus (GEO) under XXXX). Raw cyanine intensity data were processed and statistically analysed in limma [START_REF] Smyth | Linear models and empirical bayes methods for assessing differential expression in microarray experiments[END_REF]). Red and green background intensities were corrected by the o e p -method, using an offset of 50 [START_REF] Ritchie | A comparison of background correction methods for two-colour microarrays[END_REF]. A within-and betweena a o alisatio loessa d A ua tile -method, respectively) was performed on the RG data sets. Due to the different genetic background of the strains used in the two parallel experiments, all data was normalized and analysed within the experiment. The quality of the two normalized MA data sets was assessed using arrayQualityMetrics [START_REF] Kauffmann | arrayQualityMetrics-a bioconductor package for quality assessment of microarray data[END_REF]). In the linear modelling of the data, intra-spot correlations were incorporated (Smyth and Altman 2013). Significantly differentially expressed genes were identified by an empirical Bayes approach with cut-offs of the Benjamin-Hochberg adjusted p-values and log2FC at 0.05 and 0.585, respectively[e.g. used in (De Roissart, Wybouw et al. 2016)]. Only 21 and no genes exhibited significant differential transcription between the core and edge populations of the NMP and MIX experiment, respectively. Moreover, within the NMP treatment only one gene (tetur06g00230) coding for a hypothetical protein displayed an absolute fold change of higher than 2 in the edge vs core population. This is in sharp contrast to other studies where mites have exhibited strong transcriptomic responses upon genetic adaptation (Wybouw et