Giacomo Zilio
email: giacomo.zilio@umontpellier.fr.

Louise S Nørgaard

Claire Gougat-Barbera

Matthew D Hall

Emanuel A Fronhofer

Oliver Kaltz

Travelling with a parasite: the evolution of resistance and dispersal syndrome during experimental range expansion

Keywords: Dispersal syndromes, eco-evolutionary dynamics, experimental evolution, host-parasite interactions, parasitism, range expansions, resistance

Spatial dynamics of range-shifting species can be deeply affected by biotic interactions. One ubiquitous type of biotic interaction involves parasites. These can affect nearly all biological systems and impose major selective pressure on the host, leading to rapid evolutionary responses. Despite the potentially large impact of parasites, their role on host dispersal and subsequent range expansions remains mostly unexplored.

Therefore, we investigated whether parasites affect and alter host evolution during experimental range expansions. Using microbial model organisms spreading in microcosm landscapes, we found multi-trait evolution and rapid evolutionary shifts in dispersal syndromes due to spatial dynamics and parasitism. As predicted by theory, hosts that had evolved in the absence of parasites changed their movement pattern and increased dispersal at the range margins. The presence of parasites during the range expansion reshaped host phenotypic divergence between front and core, with hosts exhibiting overall reduced dispersal but increased resistance in the front. We suggest that the evolved differences in resistance and other host traits may be associated with a trade-off between dispersal and foraging efficiency.

Our work shows that the eco-evolutionary interactions between host and parasite during range expansions can shift the contenders to novel evolutionary trajectories and result in unexpected evolutionary outcomes.

Understanding and finding general patterns to these complex dynamics is of critical relevance for conservation and disease management.

Introduction

Species are on the move, and their ability to shift ranges and track climate changes have become crucial to escape the risk of extinction [START_REF] Thomas | Extinction risk from climate change[END_REF] and respond to the climate crisis and the anthropogenic alteration of the environment [START_REF] Hill | Climate and habitat availability determine 20th century changes in a butterfly's range margin[END_REF][START_REF] Parmesan | A globally coherent fingerprint of climate change impacts across natural systems[END_REF]. Additionally, biological range expansions can affect biodiversity, community structure and ecosystem functioning [START_REF] Hastings | The spatial spread of invasions: new developments in theory and evidence[END_REF][START_REF] Pecl | Biodiversity redistribution under climate change: Impacts on ecosystems and human well-being[END_REF]. Despite the importance of such dynamics, it is rarely considered how biotic interactions such as parasitism may affect the final outcome (Kubish et al. 2014). The expanding species may carry or encounter a parasite along the way and establish a new (co)evolution history, which may determine or alter the spatial spread. This could potentially lead to disease outbreaks and unpredictable results, with great concern for human health and agriculture [START_REF] Poulin | Invasion ecology meets parasitology: Advances and challenges[END_REF].

The evolution of dispersal [START_REF] Bowler | Causes and consequences of animal dispersal strategies: relating individual behaviour to spatial dynamics[END_REF][START_REF] Clobert | Dispersal Ecology and Evolution[END_REF][START_REF] Ronce | How Does It Feel to Be Like a Rolling Stone? Ten Questions About Dispersal Evolution[END_REF] and dispersal syndromes, the covariation between dispersal and other dispersal-related traits [START_REF] Clobert | Informed dispersal, heterogeneity in animal dispersal syndromes and the dynamics of spatially structured populations[END_REF][START_REF] Clobert | Dispersal Ecology and Evolution[END_REF][START_REF] Cote | Evolution of dispersal strategies and dispersal syndromes in fragmented landscapes[END_REF][START_REF] Stevens | A comparative analysis of dispersal syndromes in terrestrial and semiterrestrial animals[END_REF], is a well-established phenomenon highly supported by theoretical and growing empirical evidence (reviewed in [START_REF] Kubisch | Where am I and why? Synthesizing range biology and the eco-evolutionary dynamics of dispersal[END_REF]. It has the potential to accelerate range expansions between core and front populations, alter metapopulation dynamics and generate unexpected eco-evolutionary feedbacks [START_REF] Govaert | Eco-evolutionary feedbacks-Theoretical models and perspectives[END_REF][START_REF] Hanski | Eco-evolutionary dynamics in a changing world[END_REF][START_REF] Legrand | Eco-evolutionary dynamics in fragmented landscapes[END_REF]. However, how parasitism, the most common and ubiquitous form of biotic interaction, might affect the evolution of such dispersal syndromes remain largely unexplored. In fact, the increase of dispersal due to spatial selection, spatially assortative mating or kin competition has been usually considered in single species, e.g. bacteria [START_REF] Koskella | Using experimental evolution to explore natural patterns between bacterial motility and resistance to bacteriophages[END_REF][START_REF] Taylor | Selection experiments reveal trade-offs between swimming and twitching motilities in Pseudomonas aeruginosa[END_REF], protists [START_REF] Fronhofer | Eco-evolutionary feedbacks during experimental range expansions[END_REF], nematodes [START_REF] Friedenberg | Experimental evolution of dispersal in spatiotemporally variable microcosms[END_REF], several arthropods species [START_REF] Alford | Comparisons through time and space suggest rapid evolution of dispersal behaviour in an invasive species[END_REF][START_REF] Leotard | Range Expansion Drives Dispersal Evolution In An Equatorial Three-Species Symbiosis[END_REF][START_REF] Lombaert | Rapid increase in dispersal during range expansion in the invasive ladybird Harmonia axyridis[END_REF][START_REF] Ochocki | Rapid evolution of dispersal ability makes biological invasions faster and more variable[END_REF][START_REF] Petegem | Kin competition accelerates experimental range expansion in an arthropod herbivore[END_REF][START_REF] Simmons | Changes in dispersal during species' range expansions[END_REF][START_REF] Weiss-Lehman | Rapid trait evolution drives increased speed and variance in experimental range expansions[END_REF][START_REF] Williams | Rapid evolution accelerates plant population spread in fragmented experimental landscapes[END_REF], vertebrates [START_REF] Alford | Comparisons through time and space suggest rapid evolution of dispersal behaviour in an invasive species[END_REF][START_REF] Phillips | Invasion and the evolution of speed in toads[END_REF] and even in plants [START_REF] Szűcs | Rapid adaptive evolution in novel environments acts as an architect of population range expansion[END_REF]. Only a limited number of studies explored the consequences of species interactions on host dispersal evolution. These were focused on predation [START_REF] Pillai | Evolution of Dispersal in a Predator-Prey Metacommunity[END_REF][START_REF] Poethke | Predator-Induced Dispersal and the Evolution of Conditional Dispersal in Correlated Environments[END_REF], or on the epidemiological and ecological outcome of parasites altering the biological relationships once arrived the new habitat [START_REF] Dunn | Indirect effects of parasites in invasions[END_REF][START_REF] Strauss | Invading with biological weapons: the importance of diseasemediated invasions[END_REF]. Hence, we still have limited to no knowledge of the effect of parasites on the evolution of the spreading host and on their impact on the range expansion itself [START_REF] Kubisch | Where am I and why? Synthesizing range biology and the eco-evolutionary dynamics of dispersal[END_REF].

Rapid evolution can deeply affect the outcomes of range expansions [START_REF] Williams | Rapid evolution accelerates plant population spread in fragmented experimental landscapes[END_REF][START_REF] Williams | How Evolution Modifies the Variability of Range Expansion[END_REF], and parasites impose additional strong selective pressures on the host due to virulence and mortality [START_REF] Daversa | Infections on the move: how transient phases of host movement influence disease spread[END_REF][START_REF] Sheldon | Ecological immunology: costly parasite defences and trade-offs in evolutionary ecology[END_REF]. The little theory developed shows that parasitism may either promote or confine range expansion. Strong oscillation in host population densities due to parasitism promoted host dispersal and bet-hedging strategies, independently of dispersal costs (Chaianunporn and Hovestadt, 2012a). Interestingly, different parameter combinations not leading to such spatio-temporal variability and oscillations may select against dispersal (Chaianunporn and Hovestadt, 2012a). In another theoretical work, parasites (natural enemies and parasitoids in the paper) enforced the limits of range expansion of their hosts [START_REF] Hochberg | Can natural enemies enforce geographical range limits?[END_REF]. Similarly, modelling on plant systems emphasizes that long seed dispersal can be achieved when parasitism is associated with high host survival (low virulence) and lower dispersal costs [START_REF] Muller-Landau | Theoretical Perspectives on Evolution of Long-Distance Dispersal and the Example of Specialized Pests[END_REF]. Yet, this is rarely the case, and the ubiquity of parasites and pests [START_REF] Packer | Soil pathogens and spatial patterns of seedling mortality in a temperate tree[END_REF] impose high costs, limiting and constraining dispersal. Hosts likely rely on integrated genetic, physiological, functional or behavioural responses to face parasites during dispersal [START_REF] Brown | Stress and immunity at the invasion front: a comparison across cane toad (Rhinella marina) populations[END_REF]. Further, dispersal is a heritable trait itself [START_REF] Saastamoinen | Genetics of dispersal[END_REF] which allows population to respond to the potential differences in selection occurring in core or front populations and may drive phenotypic divergence with concurrent changes in other traits. Depending on the genetic architecture of host traits, this may or may not impose constraints on the constitutive evolutionary responses [START_REF] Hall | The Evolutionary Consequences of Stepwise Infection Processes[END_REF].The joint action of multiple selection pressures on dispersal, life-history, or interacting traits such as resistance, may lead to the evolution of different trait associations, and thus different syndromes. For example, in the case of physiological or life history trade-offs, dispersal selection at the range front may be impeded by parasite-mediated selection, such that the presence of parasites would effectively slow down a range expansion. Thus, how do the different selective pressures driving the emergence of dispersal syndromes and the response to parasitism interact in shaping host multi-trait evolution and range expansion? Through an experimental evolution approach, our goal was to investigate host phenotype diverge during a range expansion scenario and whether this was additionally affected by parasitism. We mimicked range expansions in laboratory microcosms, and in a common garden experiment we measured six relevant host traits using the host Paramecium caudatum and its bacterial parasite Holospora undulata. For simplicity we only tracked range core and range front populations: the front of infected or uninfected populations constantly dispersed into new microcosms, while the core always remained in place. We expected increased dispersal in the front, and this to be a main driver of multi-trait phenotypic divergence with concurrent changes in movement patterns, growth and population size. We predicted parasites to affect and modify host dispersal syndromes and phenotypic divergence by reducing dispersal in the infected populations and driving selection for increased resistance [START_REF] Koskella | Using experimental evolution to explore natural patterns between bacterial motility and resistance to bacteriophages[END_REF]. Since dispersal and resistance are predicted to be both costly [START_REF] Bonte | Costs of dispersal[END_REF][START_REF] Schmid-Hempel | Variation in immune defence as a question of evolutionary ecology[END_REF], the interactions and interplay of these two with the other traits, and the relative outcome of the spatial dynamic, might then depend on which selection acts stronger. Overall, the results are in accordance with our predictions confirming trait divergence between core and front populations. Parasitism modified the outcome of the range expansion reshaping host phenotypic divergence in the core and front. Considering how hosts rarely disperse alone, and that they can encounter and establish novel parasitic associations during their spread, understand the effect of parasitism on range expansion has relevant implications for biological control, conservation and management decisions.

Material and methods

Study system

Paramecium caudatum is a freshwater filter-feeding ciliate from the Northern Hemisphere [START_REF] Wichterman | The Biology of Paramecium[END_REF]. Nuclear dimorphism is typical of ciliates: The "germ-line" micronucleus is active during the sexual stage, while the highly polyploid "somatic" macronucleus regulates gene expression during the asexual stage, when replication occurs through mitotic division. In this experiment, clonal populations are maintained asexually (max. 1-2 population doublings per day at constant 23°C) in 50 ml Falcon TM tubes, using a sterilised lettuce medium (1g dry weight of organic lettuce per 1.5l of Volvic mineral water), supplemented ad libitum with the bacterium Serratia marcescens as a food resource (referred to as bacterised medium, hereafter; see [START_REF] Nidelet | Direct and correlated response to selection in a host-parasite system: testing for the emergence of genotype specificity[END_REF]. The gram-negative bacterium Holospora undulata is an obligate parasite, infecting the micronucleus of the P. caudatum [START_REF] Görtz | Diversity of Holospora bacteria in Paramecium[END_REF]. The infection life cycle comprises both horizontal and vertical transmission [START_REF] Fokin | Bacterial Endocytobionts of Ciliophora and Their Interactions with the Host Cell[END_REF] . Paramecium ingest infectious forms from the aquatic environment, which subsequently colonise the micronucleus and differentiate into multiplying reproductive forms; these reproductive forms are vertically transmitted to the daughter cells of mitotically dividing hosts. The infection life cycle is completed when reproductive forms differentiate into infectious forms, which are then released during host cell division or upon host death [START_REF] Nidelet | Direct and correlated response to selection in a host-parasite system: testing for the emergence of genotype specificity[END_REF]. Infection reduces cell division and survival of the Paramecium [START_REF] Restif | Condition-dependent virulence in a horizontally and vertically transmitted bacterial parasite[END_REF], as well as dispersal [START_REF] Fellous | Parasitic infection reduces dispersal of ciliate host[END_REF]. Experimental evolution of resistance to this parasite was demonstrated in previous long-term experiments [START_REF] Lohse | Experimental evolution of resistance in Paramecium caudatum against the bacterial parasite Holospora undulata[END_REF] and can come at reproductive costs [START_REF] Duncan | Reverse Evolution: Selection Against Costly Resistance in Disease-Free Microcosm Populations of Paramecium Caudatum[END_REF].

Long-term range expansion experiment

Dispersal arenas

Similar to Fronhofer & Altermatt (2015), we used two-patch arenas for this selection experiment (Supplementary Information, Fig. S1). The arenas were built from two 14 mL plastic tubes ("core patch" and "front patch") interconnected by 5-cm silicon tubing (0.6 mm inner diameter) serving as a corridor through which the Paramecium can actively swim and disperse. We define dispersal as the displacement of P. caudatum from the core patch to the front patch. In the long-term experiment, short episodes of dispersal (3h) alternated with periods of population growth and maintenance (1 week). Prior to each dispersal episode, the core patch was filled with 8 mL of Paramecium culture topped up with 5 mL of bacterised medium, whereas the front patch only contained 13 mL of bacterised medium (for details of the protocol, see the Supplementary Information). After the removal of a clamp that blocked the connection between the two tubes, Paramecium could freely disperse to the front patch or to stay in the core. After three hours, we blocked the corridor and estimated the cell density in the core and front patch, by sampling up to 1mL from each tube and counting the number of individuals under a dissecting microscope. The dispersal rate is thus the number of dispersers divided by the total number of individuals in the arena, divided by 3 hours.

Range expansion treatment

Two selection treatments were imposed. In the front selection treatment, only Paramecium that had dispersed into the front patch were maintained and allowed to grow for 1 week until the next episode of dispersal. Conversely, in the core selection treatment, only the non-dispersing Paramecium were maintained and allowed to regrow. These contrasting selection protocols were continued for a total of 26 cycles. The front selection treatment mimics the leading front of range expansion or a biological invasion, with populations continuously dispersing into a new microcosm. Populations from the core selection treatment stay in place and continuously lose emigrants. Each new growth cycle was started by placing on average 200 Paramecia from front and core selection treatments in 20 mL of fresh bacterised medium, carrying capacity was then reached within the following 3-4 days. The experiment was conducted with a single host line (63D). This line had undergone three years of parasite-free core selection prior to the present experiment; initially started from a mix of strains, it has become fixed for a single haplotype (O. Kaltz, unpublished data).

Parasite treatment

Core and front treatments were established for both infected and uninfected populations. Starting from a 63D laboratory culture, infected and uninfected selection lines were established. The parasites were taken from an ongoing experiment (Nørgaard et al. 2020, in prep.) that had already been imposing core and front selection on infected 63D populations for about 8 months (30 cycles). Using standard protocols (e.g., [START_REF] Duncan | Reverse Evolution: Selection Against Costly Resistance in Disease-Free Microcosm Populations of Paramecium Caudatum[END_REF] and Supplementary Information), we extracted infectious forms of the parasite from 5 core selection lines and from 5 front selection lines, which were then used to inoculate our new, naive 63D hosts.

In other words, we continued core and front selection treatments for the parasite, but replaced the previous hosts by new unselected hosts. In addition to these 10 infected selection lines, we established 3 uninfected front-selection lines and 3 uninfected core-selection lines as controls.

Adaptation assay

At the end of the long-term experiment, phenotypic trait assays for Paramecium from all 16 selection lines were performed under common-garden conditions. Parasite evolution will be analysed elsewhere (Nørgaard et al. 2020, in prep.).

Singleton isolation protocol

Using a micropipette, we arbitrarily picked 4 uninfected paramecia from each selection line and placed them individually in single 1.5 mL Eppendorf tubes filled with bacterised medium, where they were allowed to grow for 2 weeks until small monoclonal lines had established (c. 7-8 asexual generations). Each monoclonal line was then split into three technical replicates and grown for a second common-garden period of 10 days in 50-mL Falcon tubes to obtain mass cultures for the phenotype assays (16 selection lines x 4 monoclonal lines x 3 technical replicates = 192 replicates, Supplementary Information Fig. S3).

Phenotypic trait assays

After the second common-garden period and relaxed selection (absence of parasite, no dispersal treatment), technical replicates from 60 monoclonal lines were available for specifically designed tests measuring 6 phenotypic traits (Supplementary Information Fig. S3).

Resistance

To measure resistance, the Paramecium were confronted with parasites from core-selection and frontselection lines. We prepared the inocula by mixing the 5 infected core and 5 infected front selection lines, and extracting infectious forms from the two mixes (for details of the extraction protocol, see Supplementary Information 1). For inoculation, c. 5000 Paramecia were placed in a volume of 25 mL in a 50-mL tube, to which we added 4.5 x 10 5 infectious forms (core-parasite or front-parasite inoculum). In this way, we set up 4-8 inoculated tubes per host selection line, balanced between the two parasite inocula (16 selection lines x 2-4 monoclonal lines x 2 technical replicates = 120 inoculated tubes). Four days post-inoculation, we fixed c. 20 individuals from each inoculated replicate with lacto-aceto-orcein [START_REF] Görtz | Diversity of Holospora bacteria in Paramecium[END_REF]) and determined the absence or presence of infection using a phase-contrast microscope (1000x magnification). We define resistance as the proportion of uninfected individuals in the sample. Preliminary analysis showed that Paramecium from the four different selection treatments did not differ in their resistance to the mixes of front or core parasites (F 3,12 = 0.44, n.s.); we therefore combined the two inoculum sources into a single "infected" category for the main analysis.

Dispersal rate

Dispersal was measured in linear 3-patch arenas (50 mL Falcon tubes; Fig. S2), where the Paramecium dispersed from the middle tube into the two outer tubes (see Supplementary Information 1 for detailed protocol). This arena configuration allowed us to use bigger volumes of culture and to obtain higher numbers of dispersers. Connections were opened for 3 h, dispersal rates were then estimated by counting the Paramecium in samples from the central tube (500 µl) and from the combined two outer tubes (3 mL). We employed technical replicates that had not been used for the resistance assay, and were kept in 30 mL of bacterised medium for several days prior to the dispersal test (1 replicate per monoclonal line = 60 dispersal tests = 2-4 tests per host selection line).

. CC-BY-NC 4.0 International license perpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in

The copyright holder for this this version posted January 29, 2020. ; https://doi.org/10.1101/2020.01.29.924498 doi: bioRxiv preprint 8

Population growth rate and equilibrium density

For the population growth assay, we placed groups of 5 arbitrarily picked Paramecium in 15-mL tubes filled with 10 mL of bacterised medium. Over 9 days, we tracked densities in 24-h intervals, estimated from number of individuals present in 200-µL samples. We set up 6-12 tubes per host selection line (3 tubes per monoclonal line), with a total of 180 tubes tested. For each tube, estimates of intrinsic population growth rate (r 0 ) were obtained by fitting a Beverton-Holt population growth model to each density time series, using a Bayesian approach [START_REF] Rosenbaum | Estimating Parameters From Multiple Time Series of Population Dynamics Using Bayesian Inference[END_REF]. For certain tubes we obtained unsatisfactory fits of carrying capacity; we therefore decided to use the mean density over the second half of the assay (day 5-9) as a proxy for equilibrium density. 19 tubes failed to produce a coherent growth pattern and remained at very low density; it was not possible to fit our population growth model to these data, and the tubes were therefore excluded from analysis.

Swimming speed and tortuosity

At the end of the above population growth assay, we analysed swimming behaviour, using an established pipeline of computer vision and automated video analysis to collect this data [START_REF] Altermatt | Big answers from small worlds: a user's guide for protist microcosms as a model system in ecology and evolution[END_REF][START_REF] Pennekamp | BEMOVI, software for extracting behavior and morphology from videos, illustrated with analyses of microbes[END_REF]. From a given tube, one sample of 119 µL was imaged under a Perfex Pro 10 stereomicroscope, using a Perfex SC38800 camera (15 frames per second; duration: 10 s; total magnification: 10x). Videos were analysed using the bemovi R-package [START_REF] Pennekamp | BEMOVI, software for extracting behavior and morphology from videos, illustrated with analyses of microbes[END_REF] for settings see script in supplementary materials), which provided individual-based data on swimming speed and the tortuosity of swimming trajectories (= standard deviation of the turning angle distribution). Swimming speed and tortuosity were averaged over all individuals in a sample prior to analysis. A total of 60 samples (1 per monoclonal line) was used for analysis, giving 2-4 observations per host selection line.

Statistical analysis

Our main focus was the analysis of trait associations. To this end, we constructed a data matrix with the measurements of the 6 traits for 60 monoclonal lines (for resistance, the mean over the two technical replicates was calculated). To impute 10 missing observations in population growth rate and equilibrium density we used the "missMDA" package version 1.16 [START_REF] Josse | missMDA: A Package for Handling Missing Values in Multivariate Data Analysis[END_REF]. Using the approach of [START_REF] Stoks | Resurrecting complexity: the interplay of plasticity and rapid evolution in the multiple trait response to strong changes in predation pressure in the water flea Daphnia magna[END_REF], trait distributions were normalised by picking the best transformation ("bestNormalize" package version 1.4.2 of [START_REF] Peterson | Ordered quantile normalization: a semiparametric transformation built for the cross-validation era[END_REF], centred to their mean and scaled by their standard deviation. We then performed a multivariate analysis of variance (MANOVA), with range expansion treatment (core vs front selection lines), parasite treatment (infected vs uninfected control selection lines) and their interaction considered as fixed factors. Selection line identity was included as random factor. From this same data set, we also performed univariate analyses for each of the 6 traits. In a second step, to better understand the multivariate results, we performed a principal component analysis (PCA), based on the means per selection line and thus a data matrix of 16 x 6 = 96 observations. To complement graphical inspection of 9 results, we used the first and second components (PC 1, PC 2) for ad-hoc comparisons between treatments.

All statistical analyses were performed with R v 3.6.2 (R CoreTeam 2017).

Results

Univariate analyses and MANOVA

We found signatures of selection history in the observed phenotypic trait variation and covariation (Table 1).

In all univariate analyses, there were significant effects of range expansion treatment (4 traits), parasitism treatment (3 traits) or their interaction (3 traits). Except for the two swimming traits, we mostly detected simultaneous signals of the two selection treatments (Table 1, Fig. 2). This multi-trait divergence was confirmed by the MANOVA (Table 1). Looking at all traits combined, this analysis revealed significant effects of the parasitism treatment (p = 0.002) as well as significant effects of the range expansion treatment (p < 0.001) and their interaction (p = 0.046). This revealed that both selection treatments acted jointly to produce a strong phenotypic differentiation in the Paramecium.

Table 1 Results of the MANOVA and ANOVAs for the effect of range expansion and parasitism treatment and of their interaction on the six host traits analysed. The significant p-values are highlighted in bold.

Range expansion

Principle component analysis (PCA)

By means of PCA, the (co)variation in multidimensional trait space can be projected onto two main axes (PC 1 and PC 2). This allows us to describe the divergence of treatments in 2-dimensional space and to identify the individual traits that contribute most to this divergence (Fig. 1). First, along the PC 1 axis, we observed a clear effect of the range expansion treatment (front vs core selection), as illustrated by the separation of the respective clouds of points (red vs blue). This front/core separation is complete for the parasite-free selection lines (lower two clouds in Fig. 1, F 1,4 = 17.70, p= 0.013), and for selection lines evolving in the presence of the parasite (upper two clouds, F 1,4 = 6.77, p= 0.031). The direction and length of the different arrows in Fig. 1 show that the observed patterns were mainly driven by differences in swimming behaviour (PC 1 highest two loadings: speed -0.5295, tortuosity +0.5221). As found in the univariate analyses (Table 1), Paramecium from the front treatment had a 12% lower swimming speed than Paramecium from the core (Fig. 2E), instead swam in a more non-linear fashion (23% higher tortuosity; Fig. 2F). Second, along the PC 2 axis, phenotypic differentiation was driven by the parasite selection treatment. This separation between parasitised and parasite-free selection lines was complete for the front-(upper vs lower red clouds, Fig. 1, F 1,4 = 7.32, p= 0.035) but less pronounced for the core-selection treatment (upper vs lower blue points, Fig. 1, F 1,4 = 5.12, p= 0.064). Direction and length of their arrows indicate that resistance and dispersal are the key traits responsible for this multivariate response (PC 2 highest two loadings: dispersal -0.6605 resistance +0.6118).

Specifically, the univariate analyses (Table 1) show that exposure to parasites was linked with an increase in resistance for front selection lines (+ 10%, Fig. 2B), and a general decrease in dispersal (-50%, Fig. 2A) In summary, the treatment combining new front selection conditions with new exposure to parasites produced the strongest phenotypic divergence relative to the negative control treatment (no parasite, continued core selection). Multiple traits contribute to this divergence, with strongest signals from parasite resistance, dispersal and movement-related traits. The details of the PCA are provided in the Supplementary material Table S1 andS2.

Discussion

Natural populations are currently forced to climate and human-induced range-shifting during which they encounter new selective pressures as parasitism. We observed fundamental phenotypic divergence in the core and front populations during range expansions, suggesting an integrated evolutionary response of the host to spatial dynamics. Parasitism altered and reshaped such divergence, adding selection for resistance and moving host populations towards different phenotypic patterns. All host traits showed evolved differences to range expansion or parasitism, confirming the emergence of dispersal syndromes and rapid evolutionary responses to strong abiotic and biotic selective pressures.

Understanding how ecological dynamics such as range expansion affect the evolution and the trajectory of phenotypic traits is of major interest and relevance for conservation, control infectious disease and epidemic outbreaks, and to establish proper species invasion management. Range expansions are mainly driven by dispersal, a complex and multidimensional trait [START_REF] Saastamoinen | Genetics of dispersal[END_REF], determined by several characters; physiological, morphological and behavioural [START_REF] Clobert | Dispersal Ecology and Evolution[END_REF]. In line with theoretical and empirical results on range expansion dynamics, we found evidence for dispersal syndromes with increased dispersal in the range front compared to the ancestral populations, and modification in movement pattern [START_REF] Cote | Evolution of dispersal strategies and dispersal syndromes in fragmented landscapes[END_REF][START_REF] Fronhofer | Eco-evolutionary feedbacks during experimental range expansions[END_REF][START_REF] Kubisch | Where am I and why? Synthesizing range biology and the eco-evolutionary dynamics of dispersal[END_REF][START_REF] Phillips | Invasion and the evolution of speed in toads[END_REF][START_REF] Stevens | A comparative analysis of dispersal syndromes in terrestrial and semiterrestrial animals[END_REF][START_REF] Williams | How Evolution Modifies the Variability of Range Expansion[END_REF]. Counterintuitively, the high dispersers evolved lower swimming speed. This a repeatable and known result observed in other experiments with this strain of P. caudatum (Kaltz et al. individuals to change directions while swimming, which could be considered as an exploratory behaviour that increase the probability to find a dispersal corridor.

The dispersal syndromes and phenotypic divergence observed in the host primarily emerged from our experimental design of simulated range expansion. However, the arrival of the parasite imposed an additional new strong selection on the host, involving resistance and dispersal as key traits. Resistance is major trait selected under parasitism, and it has often been observed to evolve in this, and during other hostparasite evolutionary experiments [START_REF] Boots | Trade-Offs with Resistance to a Granulosis Virus in the Indian Meal Moth, Examined by a Laboratory Evolution Experiment[END_REF][START_REF] Brockhurst | The impact of migration from parasite-free patches on antagonistic host-parasite coevolution[END_REF][START_REF] Haag | Parasite-mediated selection in experimental metapopulations of Daphnia magna[END_REF][START_REF] Lennon | Is there a cost of virus resistance in marine cyanobacteria?[END_REF][START_REF] Lohse | Experimental evolution of resistance in Paramecium caudatum against the bacterial parasite Holospora undulata[END_REF]. In our work, increased resistance only evolved in the front under parasitism, even though we would have expected a similar higher level of resistance for core facing parasites compared to the ancestral stage. Additionally, we observed an overall reduced dispersal, but selection seemed to be strong enough to maintain higher level of dispersal in the front population under parasitism compared to the core, confirming parasitism to differently interact with front and core during the range expansion dynamic.

Results for a correlation of resistance and dispersal are limited and unclear [START_REF] Taylor | Bacterial motility confers fitness advantage in the presence of phages[END_REF]. Our results may suggest that a trade-off between the costly dispersal and resistance [START_REF] Bonte | Costs of dispersal[END_REF][START_REF] Schmid-Hempel | Variation in immune defence as a question of evolutionary ecology[END_REF] may emerge in the front of a range expansion (Fig. 3). The response of Pseudomonas aeruginosa to phages similarly suggested the presence of a negative trade-off between the two traits [START_REF] Whitchurch | Escherichia coli contains a set of genes homologous to those involved in protein secretion, DNA uptake and the assembly of type-4 fimbriae in other bacteria[END_REF]. In contrast, the experimental evolution of the bacterial host Pseudomonas syringae with lytic phage showed no such clear link [START_REF] Koskella | Using experimental evolution to explore natural patterns between bacterial motility and resistance to bacteriophages[END_REF]. As systemic immunological in response to infection or recovery from infections are not known in the P. caudatum-H. undulata system, another possible explanation is that resistance is associated to a negative trade-off between dispersal and foraging. P. caudatum is a filter-feeder ciliate and get infected with the bacterial parasite H. undulata by ingesting spores. If the paramecia in the front have higher dispersal and lower foraging success as we found, they might have a reduced probability of ingesting the free infectious forms in the aquatic environment.

Despite this remain speculation, the trade-off between dispersal and foraging efficiency further illustrates and well describes the evolved differences in 𝑁 ̅ and r 0 . Interestingly, our findings (high r 0 and low 𝑁 ̅ in the core vs low r 0 and high 𝑁 ̅ in the front) match a microcosm range expansion experiment with the protist Tetrahymena pyriformis [START_REF] Fronhofer | Eco-evolutionary feedbacks during experimental range expansions[END_REF]. Using a consumer-resource framework, an ecoevolutionary feedback loop based on a trade-off between dispersal and foraging was disclosed; higher dispersal and lower foraging success in the front, lower dispersal and higher foraging success in the core.

Parasitism, by impeding dispersal, may affect and modify the pace of a range expansion. This may represent either the case of an expanding species spreading with a parasite already present in the population, or the establishment of a novel host-parasite association during the spatial shift. Considering how dispersal . CC-BY-NC 4.0 International license perpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for this this version posted January 29, 2020. ; https://doi.org/10.1101/2020.01.29.924498 doi: bioRxiv preprint 14 syndromes involve the divergence and constrain on many phenotypic traits as we detected, lower dispersal may not necessarily slow down range expansion. For example, depending on the spatio-temporal habitat variability [START_REF] Poethke | Local Extinction and the Evolution of Dispersal Rates: Causes and Correlations[END_REF] or habitat niche width [START_REF] Chaianunporn | Concurrent evolution of random dispersal and habitat niche width in host-parasitoid systems[END_REF], many other traits, including life-history, may favour the colonisation and establishment of a new patch and accelerate the process of range expansion [START_REF] Burton | Trade-offs and the evolution of life-histories during range expansion: Evolution during range expansion[END_REF]. The selective pressure due to range expansion and parasitism likely modifies the genetic correlation or the genetic architecture underlying host phenotypic traits. These pressures and their interaction will produce diverse patterns of divergence, leading to changes in the whole host phenotype with unexpected traits modification or the emergence of trade-offs. The phenotype will then move into the evolutionary space and respond to new ecological selective forces, which could create ecological feedbacks direct on other traits and finally influence back evolution. The phenotypic integration is therefore essential to shed light on how complex traits like dispersal and its relationships with other traits respond to eco-evolutionary pressure in natural populations.

Here, we showed rapid phenotypic divergence of core and front populations, which was also affected and mediated by parasitism. Since these dynamics are increasing in frequency and since parasites can i) alter the success and speed of range expansions, ii) highly affect the evolution of their hosts during these spatial dynamics leading to unexpected evolutionary trajectories and eco-evolutionary feedbacks, but also iii) impact community composition, more effort towards obtaining a better comprehension and identify general patterns to this important phenomenon.

Figure 1

 1 Figure 1 PCA of the six host traits considering the first two principal component axes, PC 1 and PC 2. The length of the thin arrows inside the graph represents the loading values; the longer the arrow the higher the correlation with the PC axis and the variance explained by that trait. The points are the phenotypic values of the selection lines in multivariate space: core (blue) and front (red) of a range expansion in the presence (full dots) or absence (empty dots) of parasites. The blue empty dots with dotted ellipsis represent the ancestral host populations. The ellipses represent the 68% containment probability region of the multivariate space for the four evolutionary treatments. The black arrows on the x and y axis show the evolutionary direction of the host phenotype due to the selection imposed by range expansion and parasitism.

Figure 2 Figure 3

 23 Figure 2 Univariate response for the six traits, the grey area can be considered as the original condition of the trait. The small dots represent the mean values for the replicate lines, the big dots with the bars are the overall mean of the trait and the associated standard error. (A) Dispersal rate as the number of individuals dispersing per hour, (B) resistance as the proportion of uninfected individuals, (C) speed as µm per second, (D) tortuosity as the standard deviation of the turning angle distribution, (E) r 0 as the intrinsic rate of increase per hour, and (F) 𝑁 ̅ as the individuals per mL.

  unpublished). Moreover, speed reduction was associated with increased tortuosity, i.e. the ability of the -NC 4.0 International license perpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint inThe copyright holder for this this version posted January 29, 2020. ; https://doi.org/10.1101/2020.01.29.924498 doi: bioRxiv preprint 13

Table S1

 S1 Eigenvalues, variance explained.

		PC1	PC2	PC3	PC4	PC5	PC6
	Eigenvalues	2.4320	1.3362	1.2316	0.4925	0.2986	0.2087
	Proportion variance explained	0.4053	0.2227	0.2052	0.0820	0.0497	0.0347
	Cumulative proportion	0.4053	0.6280	0.8333	0.9154	0.9652	1
	Table S2 Loading values of the PCA.
		PC1	PC2	PC3	PC4	PC5	PC6
	Dispersal	+0.0889 -0.6605 -0.4927 +0.1638 +0.5338 +0.0324
	Resistance	+0.2123 +0.6111 -0.4568 -0.3837 +0.4015 +0.2531
	r 0	-0.4113 +0.1531 -0.5523 +0.4636 -0.4148 +0.3393
	𝑁 ̅	+0.4741 +0.2200 +0.3007 +0.7005 +0.2377 +0.2984
	Speed	-0.5295 +0.3244 +0.0500 +0.3043 +0.4754 -0.5414
	Tortuosity	+0.5221 +0.1137	-0.3877 +0.1660 -0.3148 -0.6613
			+

Acknowledgements

This work was funded by the Swiss National Science Foundation (grant no. P2NEP3_184489) to GZ and by the 2019 Godfrey Hewitt Mobility Award granted to LN by ESEB. This is publication ISEM-YYYY-XXX of the Institut des Sciences de l'Evolution -Montpellier.

Author contributions

OK, LN and GZ conceived the study. GZ, LN, CGB and OK performed the experimental work. GZ, LN, OK and EAF performed the statistical analysis. All authors interpreted the results. GZ and OK wrote the first draft of the manuscript and all authors commented on the final version.

. CC-BY-NC 4.0 International license perpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for this this version posted January 29, 2020. ; https://doi.org/10.1101/2020.01.29.924498 doi: bioRxiv preprint

Supplementary Information

Extraction protocol

Inocula were prepared concentrating by centrifugation infected paramecia (35000 RPM for 20 minutes) in Falcon tubes filled with 15 mL of medium. After removing the supernatant, the concentrated individuals were transferred into 1.5 mL Eppendorf tubes containing 1 mm glass beads, and they were vortexed and crushed using a Qiagen TissuLyser (1.45 minutes at 30 oscillation frequency). The released infectious forms were then counted using a hemocytometer at 200 x magnification under a microscope (Leica DM LB2), and their concentration was adjusted with sterile water.

Dispersal arenas

We first filled the two-patch system with 9.5 mL of fresh growth medium and we then closed the corridor with a clamp. Secondly, we filled to 12 mL one of the two tubes, the core patch, using the paramecia and medium from an experimental selection line. The second tube, the front patch, was filled to 12 mL with fresh growth medium, and thus resulted empty at this stage. Thirdly, we removed the clamp and opened the corridor allowing the paramecia to actively disperse and swim from core to front patch or to stay in the core.

After three hours, we closed the corridor tube and we estimated the population density by taking a 200 µL sample from core and front patch and counting the number of individuals under a dissecting microscope. The arenas were composed from two patches (core and front) interconnected by 5-cm silicon tubing serving as a corridor. The paramecia could displace from core to front patch through active dispersal once we removed the clamp blocking the corridor. In the long-term experiment, we alternated short episodes of dispersal (3h) with periods of population growth and maintenance (1 week = 1 cycle).

Dispersal rate

We followed a similar procedure to what was previously described for thedispersal arenas. In short, we filled the 3-patch systems with 40 mL of fresh growth medium, we closed the two corridors with clamps, and then added the entire 25 mL of the final tube to the central compartment. Left and right tubes were filled with fresh growth medium to have the same amount of volume of the central tube, then the corridors were opened. We allowed the paramecia to actively disperse in both directions for three hours before closing the dispersal corridors again. We sampled 500 µL from the core patch, 3 mL from both front patches and we estimated population densities. 500µl) and from the combined two outer tubes (3 mL). We did not control for the density of Paramecium placed in the middle tube of the dispersal arena. However, preliminary analysis showed no significant effect of density on dispersal rate (F1,54= 0.69, n.s.), and the covariate was therefore omitted from further analyses.

Figure S3

Experimental design of the adaptation assays (16 x 4 x 3 tubes). Some technical replicates were lost during the manipulation, and the final number was 180. This corresponded to the total of 60 monoclonal lines that were used in the analysis.


Julien Allali
email: allali@univ-mlv.fr

Marie-France Sagot
email: sagot@inria.fr

A Multiple Graph Layers Model with Application to RNA Secondary Structures Comparison

Keywords: graph layers, graph comparison, edit distance, RNA, secondary structure

We introduce a new data structure, called MiGaL for "Multiple Graph Layers", that is composed of various graphs linked together by relations of abstraction/refinement. The new structure is useful for representing information that can be described at different levels of abstraction, each level corresponding to a graph. We then propose an algorithm for comparing two MiGaLs. The algorithm performs a step-bystep comparison starting with the most "abstract" level. The result of the comparison at a given step is communicated to the next step using a special colouring scheme. MiGaLs represent a very natural model for comparing RNA secondary structures that may be seen at different levels of detail, going from the sequence of nucleotides, single or paired with another to participate in a helix, to the network of multiple loops that is believed to represent the most conserved part of RNAs having similar function. We therefore show how to use MiGaLs to very efficiently compare two RNAs of any size at different levels of detail simultaneously.

Introduction

We introduce in this paper a new data structure, called MiGaL for "Multiple Graph Layers", that is composed of various graphs linked together by relations of abstraction/refinement. The new structure is useful for representing information that can be described at different levels of abstraction, each level corresponding to a graph. Similar structures have already been used for modelling a spatial environment [4,[START_REF] González | Multihierarchical graph search[END_REF] or plant architectures [START_REF] Godin | A multiscale model of plant topological structures[END_REF][START_REF] Ferraro | An edit distance between quotiented trees[END_REF]. It could also be used for analysing web pages or source codes. MiGaL is a more general structure in that it can model graphs besides trees, and edges can be added between the different layers.

After giving a formal presentation of the MiGaL structure, we propose an algorithm for comparing two MiGaLs. The algorithm performs a step-by-step comparison starting with the most "abstract" level. The result of the comparison at a given step is communicated to the next step using a special colouring scheme.

We then present an application of MiGaL to the analysis of RNA secondary structures. RNAs have different functions in a cell, these being strongly related to the spatial fold adopted by the molecule. It is therefore meaningful to compare such folds in order to infer or understand function. For RNAs, the comparison is usually done by considering the secondary structure which, in the case of RNAs, provides a planar view of the fold. A secondary structure may be seen at different levels of detail, going from the sequence of nucleotides, single or paired with another to participate in a helix, to the network of multiple loops that is believed to correspond the most conserved part of RNAs having similar function.

MiGaLs represent therefore a natural way of modelling such secondary structures. To perform the comparison at each different level, an edit distance algorithm is applied to trees that are rooted and ordered (reflecting the orientation of the RNA molecule). We use for this the algorithm particularly adapted to RNAs that was introduced in a recent work by Allali et al. [START_REF] Allali | A new distance for high level rna secondary structure comparison[END_REF]. In this paper, we show how to optimise it using the colouring scheme of the algorithm for comparing two MiGaLs introduced earlier in the paper in order to compare RNAs at different levels of detail simultaneously. We show that our method leads to quite satisfying results when applied to the analysis of RNA secondary structures. Furthermore, we show that this new algorithm for comparing RNAs allows to address the scattering effect problem mentioned in the literature [START_REF] Allali | A new distance for high level rna secondary structure comparison[END_REF]. Finally, the divide strategy used as the different levels are considered in turn allows also to very efficiently compare RNAs of any size.

Multiple Graph Layers

We present in this section the MiGaL data structure. After a formal definition, we provide an algorithm to compare such structures.

Definition

From now on, we adopt the following notations. Given a set S, we denote by |S| its cardinality. Let G(V, E) be a graph with vertex set

V (G) = {v 1 , . . . , v |V (G)| } and edge set E(G) = {e 1 , . . . , e |E(G)| }.
We start by defining the new data structure. It can be described as a layered sequence of graphs, with each graph linked to its neighbours by two applications, one of which corresponds to an abstraction from the previous graph, the other a refinement leading to the next graph. The formal definition of the MultI-GrAphLayers data structure is as follows.

Definition 1 (Definition of a M iGaL). A M iGaL structure M (G, R) of size |M | is defined as a layered sequence G = G 1 , . . . , G |M| of |M | graphs and a sequence R = α 1 , . . . , α |M|-1 of |M | -1 applications called refinements. Each refinement α i is an application from V (G i ) to P(V (G i+1 )), that is, each vertex in V (G i ) has a subset of V (G i+1
) as image. The inverse application of α i , denoted by β i , is a surjection called an abstraction that maps every vertex in V (G i+1 ) to a vertex of V (G i ). In addition, M (G, R) satisfies the three conditions:

1. For each vertex v ∈ V (G i ), the subgraph induced by α i (v) is connected and not empty. 2. For each edge (u, v) ∈ E(G i ), there exists at least one vertex of α i (u) connected to a vertex of α i (v). 3. For each pair of unconnected vertices u, v ∈ V (G i ), there is no edge between a vertex of α i (u) and a vertex of α i (v). This structure is clearly useful to encode data at different levels of detail. A comparable data structure called Multihierarchical Graph has been used in [4,[START_REF] González | Multihierarchical graph search[END_REF] to represent the environment (towns, buildings, rooms, . . . ) and drive robots in space. Godin et al. in [8,[START_REF] Ferraro | An edit distance between quotiented trees[END_REF] use a similar structure, that they call quotiented trees, to model and study plants. As compared to the latter, MiGaLs can represent graphs and not just trees, and contrary to Multihierarchical Graphs, MiGaLs can have edges between nodes belonging to different layers.

Top-Down comparison algorithm

In this section, we present an algorithm to compare two MiGaLs. The algorithm performs a top-down traversal of the structure. We assume both MiGaLs have the same number of layers, and compare them layer by layer starting with the graph at the top. The result of the comparison at a given layer is transmitted to the next layer by colouring vertices and edges of the graph and using the refinement application. We then compare the two graphs of the next layer taking into account such colouring. The process continues until the last layer is reached (bottom graph).

This approach assumes the existence of an algorithm allowing to compare two graphs of a same layer. This algorithm clearly depends on the data that is modelled using a MiGaL. In what follows, we consider a black box that is able to produce an extended mapping between two graphs where by an extended mapping is meant the following: Definition 2 (Extended mapping). Given two graphs G 1 and G 2 , we define an extended mapping M between them as a set of couples of vertex sets of the two graphs:

M = {(S 1 , S 2 )/S 1 ⊂ V (G 1 ) and S 2 ⊂ V (G 2 )} such that for all (S 1 , S 2 ) ∈ M: -S 1 is a connected component of G 1 -S 2 is a connected component of G 2 -∀(S ′ 1 , S ′ 2 ) ∈ M, S 1 ∩ S ′ 1 = ∅ or S 2 ∩ S ′ 2 = ∅
We define a colour-partitioned graph G with vertex set V (G) as a graph fitted with an application C G : S ⊂ V (G) → N + which gives a colour to each vertex of S ⊂ V (G) such that subgraphs defined by vertices of the same colour are connected. For convenience, uncoloured vertices are given the colour 0. The vertices with colour 0 are therefore elements of V (G) \ S.

We now extend the definition of an extended mapping to colour-partitioned graphs as follows.

Definition 3 (Colour-constrained extended mapping). Given two colourpartitioned graphs G 1 and G 2 , a colour constrained extended mapping M between them is defined as an extended mapping which further satisfies the following conditions:

-∀(S 1 , S 2 ) ∈ M, every vertex of S 1 has the same colour c and every vertex of S 2 also has this colour c. -There is no vertex of G 1 or G 2 coloured with 0 that is involved in M. Now that we have defined a mapping on colour-partitioned graphs, we give a general algorithm for comparing two MiGaLs, M (G, R) and M ′ (G ′ , R ′ ), having the same number of layers. This number is denoted by |M | = |M ′ |. The algorithm assumes again that we have a black box B that computes a colour-constrained mapping between two colour-partitioned graphs.

The initialisation step colours the vertices of G 1 and G ′ 1 with 1. The algorithm is then divided into |M | steps. For each layer i from 1 to |M |, we compute the mapping M i between G i and G ′ i using B and (except for the last step) colour the vertices of G i+1 and G ′ i+1 such that vertices associated by the mapping have the same colour and vertices absent from M i are coloured with 0. The result of the algorithm is the set of all mappings between each pair of layers.

The pseudo-code of this algorithm is shown in Figure 2. Let b be the time complexity of algorithm B; the time required to compare two MiGaLs of size n is O(n * b).

Compare(M (G, R), M ′ (G ′ , R ′ )) 1. colour = 2 2. Set all vertices of G1 and G ′ 1 to 1. 3. for each layer i from 1 to |M | 4. M = B(Gi, G ′ i ) 5.
set the color of V (Gi+1) and V ′ (G ′ i+1 ) to 0 6.

for each couple of sets (u, v) in M 7.

for all nodes n in u 8.

set the colour of nodes in α(n) to colour 9.

for all nodes m in v 10.

set the colour of nodes in α ′ (m) to colour 11.

set colour to colour + 1 12. return M 1...|M | Fig. 2. Algorithm for the comparison of two MiGaL structures using an external algorithm B that computes a colour-constrained mapping between two colour-partitioned graphs.

It is important to notice that the main idea of this algorithm is to compute mappings from the top to the bottom layer. Each layer is compared using the mapping of the previous one without reconsidering the choices implied by this mapping.

We now present a practical application of the MiGaL structure and of the comparison algorithm to the study of RNA secondary structures. In particular, we show, in this case, how to significantly improve the complexity of algorithm B by taking advantage of the node-colouring.

RNA-MiGaL

RNAs are one of the most important molecules in the cell. They are composed by a succession of nucleotides named A,C,G and U (also referred to as bases). Inside a cell, RNAs do not keep a linear form but instead fold in space. The fold is given by the set of interactions between nucleotides. An RNA can be described at three different levels, respectively called its primary, secondary and tertiary structures. The primary structure refers to the sequence of nucleotides. The secondary structure is composed of the list of base pairs that participate in an helix (see below). The tertiary structure corresponds to all interactions (base pairings) in the RNA that is, to its 3D structure.

The function of an RNA (be it the well known ribosomal and transfert RNAs or the more recently discovered snoRNAs, microRNAs, etc.) is strongly linked to the shape adopted by the RNA in space. It is accepted that two RNAs that have the same function will have closely similar secondary structures but not necessarily similar primary structures. Considering this, it is fundamental to have efficient algorithms to compare RNAs from the point of view of their secondary structures.

Various structural elements can be distinguished in the secondary structure of an RNA: helices which correspond to consecutive base pairs, hairpin-loops which are unpaired bases at the end of an helix, internal-loops defined by the unpaired bases between two helices and called bulges when one side is empty, multi-loops which are the meeting point of at least three helices and, finally, stems that are series of helices, internal-loops and bulges.

Many approaches have been used for modelling secondary structures. Mainly, three codings has been proposed to represent RNA secondary structures: rooted ordered trees [START_REF] Zuker | RNA secondary structures and their prediction[END_REF]; arc annotated sequences [START_REF] Evans | Finding common subsequence with arcs and In the 10th Annual Symposium Combinatorial (CPM'99)[END_REF]; 2-intervals [START_REF] Vialette | Pattern matching problems over 2-interval sets[END_REF]. Fig. 3. Edit operations, edit distance and alignment on rooted ordered labelled (on nodes and edges) trees. On the first line, the first tree is edited by a substitution of a-3 to d-5. Then the node 5 is deleted and after that the node 0 is inserted. From this point, if operations have unit cost (1 for mismatch, deletion, insertion), the edit cost is 3 (the edit distance is 2 because the substitution is useless). The last two trees on the first line give an example of edge fusion (edges ce and b, the new label is constructed by concatenation) and node fusion (nodes 0 and 1, the new label is computed using mean). The second line shows how to align the first tree of the first line with the fourth. We see here that the alignment cost is four (two insertions and two deletions).

Here we focus on the tree model. We represent secondary structures by rooted ordered trees because the root correspond to the beginning and the end of the molecule and the order between the children of a node corresponds to the orientation of the sequence. In fact, we can use different tree models depending on the information we want to encode. In [START_REF] Zuker | RNA secondary structures and their prediction[END_REF], Zuker and Sankoff use trees where internal nodes code for base pairs and leaves for free nucleotides. This tree can be compacted [START_REF] Fontana | Statistics of RNA secondary structures[END_REF] into a homeomorphically irreductible tree where internal nodes correspond to helices and leaves correspond to fragments of unpaired bases.

Shapiro in [START_REF] Shapiro | An algorithm for multiple RNA secondary structures[END_REF] uses a tree where edges code for helices and nodes are labelled on {M, B, I, H, R} (multi-loop, budge, internal-loop, hairpin-loop and the root). We can also use a more abstract tree where edges code for stems and nodes for multi-loops and hairpin-loops, or just for multi-loops.

To compare rooted ordered trees, there exist essentially two methods. The first is the edit distance [START_REF] Selkow | The tree-to-tree editing problem[END_REF][START_REF] Tai | The tree-to-tree correction problem[END_REF][START_REF] Zhang | Simple fast algorithms for the editing distance between trees and related problems[END_REF] based on three edit operations. A substitution changes the label of a node. A deletion removes a node of the tree, its children are then re-attached to the node's father preserving the relative order between nodes. An insertion is the opposite operation of a deletion. If we assign a score to each of these operations, we can define the edit distance between two trees as the minimum of the score of a series of edit operations (sum of the score of each operation) that transforms the first tree into the second one. Recently, Allali et al. [START_REF] Allali | A new distance for high level rna secondary structure comparison[END_REF] extended this distance by introducing four new operations: node fusion, edge fusion and the opposite operations of node split and edge split. These operations allow to address some of the limitations of the classical edit distance when comparing RNA secondary structures using high level trees (where nodes and edges code for secondary structure elements only, not for nucleotides). The complexity of the classical edit distance is O(n 4 ) where n is the size of the trees and O((2d) l n 4 ) for the algorithm with node and edge fusions where d is the degree of the trees and l the maximum number of consecutive fusions per node.

The second method for comparing rooted ordered trees is by performing a tree alignment [START_REF] Jiang | Alignment of trees -an alternative to tree edit[END_REF][START_REF] Höchsmann | Local similarity in RNA secondary structures[END_REF]. In this case, the goal is to insert blank labelled nodes in the two trees such that they become isomorphic. The score of the alignment corresponds to the sum of the scores of the association of node labels (the blank character assumes the role of insertions and deletions in sequence alignment). The optimal score, maximal or minimal depending on the scoring scheme adopted, is then sought. A tree alignment can be expressed as an edit distance computation where insertions must precede all deletions. We do not give further details here on alignments as we only use edit distance later in the paper.

We now suggest a new modelling of RNA secondary structures based on the MiGaL data structure introduced earlier in the paper. We thus call RNA-MiGaL a MiGaL structure composed by four layers, each modelled by a rooted ordered tree. The next section is dedicated to a description of RNA-MiGaL.

The four layers definition

The four layers contained in an RNA-MiGaL correspond to the secondary structure of an RNA observed at different levels of detail. Thus, each layer is modelled by a rooted ordered labelled tree. In the bottom layer, nodes and leaves code for nucleotides while the top layer encodes the network of multi-loops of an RNA. This choice has been dictated by two assumptions. The first, already mentioned, is that structure is more important than sequence. Thus we introduce information on nucleotides at the bottom layer only which will be treated last by the comparison algorithm. The second is that the network of multi-loops can be considered as the skeleton of the secondary structure of an RNA. RNAs of the same family should thus have strongly conserved multi-loop networks. For this reason, the top layer of an RNA-MiGaL is a tree that codes for the multi-loop network. Intermediate layers correspond to the structure encoded using stems or helices.

We provide below a summary of the layers and of the meaning of a node and a leaf in each layer:

-The tree of layer 1 corresponds to the multi-loop network. The nodes encode multi-loops and the edges encode stems. In the nodes, we store the number of helices connected by the multi-loop. In the edges, we store the number of nucleotides contained in the stem. -Layer 2 consists in the structure defined by the stems. The internal nodes represent multi-loops, leaves represent hairpin-loops and edges encode stems.

In the edges, we store the number of base pairs and the number of unpaired bases contained in the stem. In the internal nodes and leaves, we store the number of unpaired bases contained in the multi-loops and hairpin-loops. -The tree of layer 3 encodes secondary structure elements: nodes encode hairpin-loops, multi-loops, internal-loops and bulges and store the number of unpaired bases of the corresponding element. The edges represent helices and store the number of base pairs of the helices. -The last tree models the RNA primary structure. The internal nodes thus represent base pairs and leaves unpaired bases. Both store the names of the corresponding bases.

RNA-MiGaL comparison

The problem now is to compare two secondary structures using RNA-MiGaLs.

To do so, we use edit distance with fusions and the algorithm described in [START_REF] Allali | A new distance for high level rna secondary structure comparison[END_REF] to compare the pairs of trees of layers 1, 2 and 3. To take into account the colour of the nodes and edges, we just have to test if two nodes or edges have the same colour to allow them to be fusioned (in a tree) or substituted (between the trees). The trees of layer 4 are compared using the classical edit distance as fusions make no sense for nucleotides. Figures 4 and5 show the result of the comparison using RNA-MiGaLs between two Group I Intron RNAs retrieved from [START_REF] Cannone | The comparative RNA web (CRW) site: an online database of comparative sequence and structure infor mation for ribosomal, intron, and other RNAs[END_REF]. The left RNA is found in Acanthamoeba griffini and the right one in Chlorella sorokiniana.

The time complexity required to compare two trees of a same layer is O((2d) l n 4 ) for the first 3 layers and O(n 4 ) for the last layer with n the size of the trees. Since we colour nodes while we progress in the comparison, and, once trees are coloured, substitutions and fusions can only be performed on nodes of same colour, we may wonder if it is not possible to compare separately the subtrees defined by each colour.

The response is positive but it requires some attention as is shown in the following example. Figure 6 presents on the top left a colour-partitioned tree, with white corresponding to colour 0. In the box, we have split each tree into subtrees according to their colour and computed the edit distance between these trees separately. The nodes associated by the mapping are linked with dashed

lines. We have then reported the result of computations

The problem is pointed by bold lines in the tree at the bottom. We have two couples of nodes associated via two different computations. Inside each computation, the relative order of the nodes is respected by the edition (we work on ordered trees). However, on the final tree, the associations do not respect the order between the nodes. We therefore have to the algorithm such as to take into account the presence of a subtree of a colour x hanging from an edge of a subtree of colour y = x. To do so, we introduce anchors during the splitting step. When we split trees into subtrees according to their colour, we add anchors to the subtrees. These anchors represent the subtrees of another colour that hang from an edge of the subtree being considered as shown in Figure 7. The anchors have the same colour as the subtrees they stand for. We then have to modify the edit distance scoring scheme such that:

deleting the anchor of colour c costs the deletion of the subtree corresponding to the anchor; -an anchor of colour c can only match with an anchor of the same colour.

Finally, we order the edit computations according to the dependence implied by the colouring scheme, that is we consider the subtree from the leaves to the root.

The consequence on the complexity is important as the time required to compare two trees of a given layer has a time (and space) complexity which depends on the size of the biggest subtree and not on the size of the whole trees anymore.

In the worst-case, the time complexity is 0((2d) 4 for each in case of RNAs, on average, the tree at each layer contains number of subtree that is proportional to the size of the tree, and the observed time complexity is close to linear except for the first layer.

Conclusion

We introduced in this paper a new data structure and an algorithm for comparing structures. Both are generic and may be applied to very different types of problems. They are particularly interesting for comparing objects that may be described at different levels of abstraction.

RNA secondary structures are one example of such objects and we therefore showed how to optimise the generic algorithm, in particular its special colouring scheme for going from one level of abstraction to the next, to compare two RNA secondary structures. The results obtained are very satisfying. In particular, the algorithm addresses the so-called scattering effect described in the literature [START_REF] Allali | A new distance for high level rna secondary structure comparison[END_REF]. The new algorithm allows to compare large sized structures in a fast way while the multiple layer approach represents a biologically very natural way of modelling and analysing with RNAs.

Fig. 1 .

 1 Fig. 1. Example of a M iGaL structure of size three defined by three graphs and two refinements.

Figure 1

 1 Figure 1 illustrates this definition. From now on, we refer to the first graph of a MiGaL as the top graph and to the last one as the bottom graph.This structure is clearly useful to encode data at different levels of detail. A comparable data structure called Multihierarchical Graph has been used in[4,[START_REF] González | Multihierarchical graph search[END_REF] to represent the environment (towns, buildings, rooms, . . . ) and drive robots in space.Godin et al. in [8,[START_REF] Ferraro | An edit distance between quotiented trees[END_REF] use a similar structure, that they call quotiented trees, to model and study plants. As compared to the latter, MiGaLs can represent graphs and not just trees, and contrary to Multihierarchical Graphs, MiGaLs can have edges between nodes belonging to different layers.

Fig. 4 .

 4 Fig. 4. Result of the comparison of two Group I Introns. On the left, the result of the comparison of the trees of layer 1; on the right, the result of the comparison of the trees of layer 2.

Fig. 5 .

 5 Fig. 5. Result of the comparison of two Group I Introns. On the top, the comparison of the trees of layer 3; on the bottom, the result of the comparison of the trees of layer 4.

Fig. 6 .

 6 Fig.6. Problem about splitting tree according to colour to optimize edition computation.

Fig. 7 .

 7 Fig. 7. Splitting the left tree into subtrees according to the colours of the nodes and addition of anchors (represented by triangles).